/*
 *  glmesh2.c:  2-D triangle mesh generation using OpenGL.
 *
 *  This program reads a file containing a 2-D triangle mesh, allows
 *  an interactive user to update the mesh with additions, deletions,
 *  and alterations, and writes a file containing the revised mesh.
 *
 *  R. Renka
 *  04/29/2009
 *
 *  Compile and link command for Linux:
 *
 *     gcc -O3 glmesh2.c -o glmesh2 -I/usr/X11R6/include -L. -lm -lglut
 *
 *  Compile and link command for Mac OS X:
 *
 *     gcc -O3 glmesh2.c -o glmesh2 --ansi -framework GLUT
 *         -framework OpenGL -framework Cocoa
 *
 *  Linux and Windows need <GL/glut.h>; OS X needs <GLUT/glut.h>.
 *
 *  Required ordering of header files (not all of which are required):
 *    windows.h, stdlib.h, glew.h, glut.h.
 * 
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#include <GL/glut.h>   /* Linux version */

struct rect {    /* Rectangle in dialog box window */
   int x0, y0;   /* Window coordinates of lower left corner */
   int w, h;     /* Width and height in window coordinates */
};

/*
 *  Static array dimension:
 */
#define LSTRING 80  /* Length of title string and dialog string */

/*
 *  Enumeration constant and arrays of strings for printing.
 */
enum Mode {DEFAULT, SELECT_CIRCLE, PAN, CONSTRUCT_R, ALTER_R,
           ALTER_VIEW_VOLUME, ALTER_DOMAIN, SWAP_EDGES,
           PLACE_TITLE, EDIT_UV};
static char *modeName[] = {"DEFAULT", "SELECT_CIRCLE", "PAN",
                           "CONSTRUCT_R", "ALTER_R", 
                           "ALTER_VIEW_VOLUME", "ALTER_DOMAIN",
                           "SWAP_EDGES", "PLACE_TITLE", "EDIT_UV"};
static char *flagName[] = {"FALSE", "TRUE"};

/*
 *  Global variables defining the triangle mesh:
 *
 *  A triangle mesh is a set of triangles in which the intersection
 *  of two triangles, if not empty, is a vertex of both or an edge 
 *  of both, and whose union is a connected region with one or more
 *  simple closed non-intersecting boundary curves.  Boundary curves
 *  are oriented so that triangles are on the left as the sequence 
 *  of boundary vertices is traversed.  The domain geometry can be 
 *  defined by designating a subset of the boundary vertices as being
 *  non-removable.  The following variables characterize the mesh.
 *
 *  nvmax = Maximum number of vertices for which storage has been
 *          allocated.  This value is doubled, and additional storage
 *          is allocated when necessary as the mesh is updated with 
 *          additional vertices.
 *
 *  nc >= 1:  Number of boundary curves.
 *  nb >= 3*nc:  Number of boundary vertices (boundary edges).
 *  nv >= nb:  Number of vertices.
 *  nn <= nv:  Number of non-removable vertices.
 *  nt = 2*nv - nb + 2*nc - 4 <= 2*nv - 5:  Number of triangles.
 *  ne = 3*nv - nb + 3*nc - 6 <= 3*nv - 6:  Number of edges.
 *
 *  vtxy = Array dimensioned nvmax by 2 containing the Cartesian
 *         coordinates of the vertices.
 *
 *  uv = Array dimensioned nvmax by 2 containing pairs of function
 *       values (u,v), such as velocity vectors, at the vertices.
 *       The values default to zeros, but may be altered.
 *
 *  ltri = Array dimensioned 2*nvmax by 9 containing a triangle list
 *         in the first nt rows:  CCW-ordered vertex indices (0 to 
 *         nv-1), neighboring triangle indices (-1 to nt-1), and
 *         opposite edge indices (0 to ne-1), where nonexistent
 *         neighboring triangles are indexed by -1.  The ordering of
 *         the neighboring triangle indices and edge indices is 
 *         defined by that of the vertices:  edge ltri[kt][i+6] is 
 *         opposite vertex ltri[kt][i] and is shared by triangles kt
 *         and ltri[kt][i+3] for i = 0, 1, and 2.
 *
 *  lvtx = Array of length nvmax containing a vertex list in the 
 *         first nv positions:  index of a triangle that contains
 *         vertex kv in position kv for kv = 0 to nv-1.  If kv indexes
 *         a boundary vertex, then triangle lvtx[kv] also contains the
 *         next vertex in the boundary curve, where boundary curves
 *         are oriented so that triangles are on the left.
 *
 *  ledg = Array of length 3*nvmax containing an edge list in the 
 *         first ne positions:  index of a triangle (the one with
 *         larger index if there are two) that contains edge ke in 
 *         position ke for ke = 0 to ne-1.  
 *
 *  lnrv = Array of length nvmax containing the indices (for vtxy) of
 *         the non-removable vertices in the first nn positions.
 *
 *  lcrv = Array of length nvmax/3 containing the index of one 
 *         boundary vertex in each boundary curve in the first nc
 *         positions.
 */
static GLint nvmax, nc, nb, nv, nn = 0, nt, ne;
static GLdouble (*vtxy)[2], (*uv)[2];
static GLint (*ltri)[9], *lvtx, *ledg, *lnrv, *lcrv;

/*
 *  Global variables associated with the currently selected region R:
 *
 *  The user may select a polygonal region R to be refined, optimized,
 *  coarsened, or removed from the mesh.  The polygon boundary is a 
 *  simple closed unoriented curve defined by a sequence of three or 
 *  more vertices (selected by mouse button presses).  The following
 *  variables are associated with R.
 *
 *  mode = Current mode (described below).  Values CONSTRUCT_R and
 *         ALTER_R are related to R. 
 *  newR = Flag with value TRUE iff Function initR needs to be called
 *         to store listf, lister, listtr, and listvr before the mesh
 *         can be altered:  set to TRUE when R is created and when R 
 *         is altered by adding or moving a vertex (Function mouse); 
 *         set to FALSE by Function initR; and tested by Functions 
 *         coarsenRv, convertRv, fillRt, moveRn, optimizeRt, 
 *         optimizeRv,  refineRe, refineRt, removeRt, and trprint.
 *  indxR = Index (for R, in the range 0 to np-1) of a vertex 
 *          currently being dragged by the mouse.
 *  np = Number of vertices in the boundary of the currently selected
 *       polygonal region R.
 *  npmax = Maximum value of np for which it is not necessary to
 *          reallocate storage for R.
 *  listf = Array of length nvmax used to store Boolean flags
 *          marking vertices as contained in R.
 *  lister = Array of length 3*nvmax containing the indices of the
 *           edges that are contained in R (first ner positions).
 *  listtr = Array of length 2*nvmax containing the indices of the
 *           triangles that are contained in R (first ntr positions).
 *  listvr = Array of length nvmax containing the indices of the
 *           (triangle mesh) vertices that are contained in R (in the
 *           first nvr positions).
 *  lste = Temporary storage array of length 3*nvmax used by Functions
 *         optimizeRt and refineRe to store a list of edges to be 
 *         swapped, and used by Functions coarsenRv, fillRt, 
 *         getCorners, dragCorner, and optimizeRv to store vertex 
 *         lists.
 *  ner = Number of edges (whose vertices are) contained in R.
 *  ntr = Number of triangles (whose vertices are) contained in R.
 *  nvr = Number of (triangle mesh) vertices contained in R.
 *  R = Array dimensioned npmax by 2 containing the Cartesian coordi-
 *      nates of the sequence of vertices defining the boundary of R.
 */
static GLboolean newR = 0;
static GLint indxR;
static GLint np = 0;
static GLint npmax = 100;
static GLboolean *listf;
static GLint *lister, *listtr, *listvr, *lste;
static GLint ner = 0, ntr = 0, nvr = 0;
static GLdouble (*R)[2];

/*
 *  Additional global variables:
 *
 *  antialias = Flag with value TRUE iff edges are to be drawn with
 *              antialiasing.
 *  aspect = Flag with value TRUE iff the viewport is required to
 *           have the same aspect ratio as the view volume.
 *  bdry_layer = Flag with value TRUE iff the smoothing function
 *               optimizeRv is to create a boundary layer by moving
 *               vertices adjacent to boundary vertices closer to
 *               the boundary.
 *  color_code = Flag with value TRUE iff triangles are to be color-
 *               coded by quality using function trqual1, as opposed
 *               to being rendered in outline mode by drawing only
 *               the edges.
 *  color_b = Color of background.
 *  color_e = Color of triangle mesh edges.
 *  color_t = Color of text strings.
 *  cr,cx,cy = Radius and center of a circle selected by Functions
 *             getCircle and motion.
 *  dlg_exit = Flag with value TRUE iff the program is to be
 *             terminated when the user selects the 'Yes' button in
 *             the dialog box.
 *  dlg_r0 = Rectangle representing 'No' button.
 *  dlg_r1 = Rectangle representing 'Yes' button.
 *  dlg_string = Text string displayed in dialog box window.
 *  dlg_yes = Flag with value TRUE iff the user selects the 'Yes'
 *            button in the dialog box.
 *  dx,dy = Width and height of view volume.
 *  fname = File name for output.
 *  ind_font = Font number for indices.
 *  line_width = Line width in pixels for drawing edges.
 *  mode = Current mode, altered by a keypress or menu selection, and
 *         tested by Function motion:  
 *         mode = SELECT_CIRCLE if the user is selecting a circle.
 *                Refer to Function moveRn.
 *         mode = PAN if the user is selecting a new center (xc,yc).
 *         mode = CONSTRUCT_R if the user is selecting a new polygonal
 *                region R.
 *         mode = ALTER_R if the user is altering polygonal region R.
 *         mode = ALTER_VIEW_VOLUME if the user is selecting a new
 *                viewing volume by dragging a rectangle.
 *         mode = ALTER_DOMAIN if the user is altering the domain by
 *                dragging non-removable vertices.
 *         mode = SWAP_EDGES if the user is selecting triangle edges
 *                in convex quadrilaterals for edge swaps.
 *         mode = PLACE_TITLE if the user is selecting a title 
 *                location.
 *  nnbv = Number of neighbors of a non-removable vertex that is 
 *         currently being dragged with mode = ALTER_DOMAIN.
 *  nref = Number of refinements (by refineRe or refineRt) that have
 *         not been reversed by Function unrefine, or rendered non-
 *         reversible by Function coarsenRv or removeRt.
 *  nrefmax = Maximum value of nref for which it is not necessary to
 *            reallocate storage for rstack.
 *  plot_box = Flag with value TRUE iff an axis-aligned bounding box 
 *             with labeled corner coordinates is to be drawn.
 *  plot_e = Flag with value TRUE iff edge indices are to be displayed
 *           (at edge centers).
 *  plot_nrv = Flag with value TRUE iff non-removable vertex indices 
 *             are to be displayed.
 *  plot_t = Flag with value TRUE iff triangle indices are to be 
 *           displayed (near barycenters).
 *  plot_v = Flag with value TRUE iff vertex indices are to be 
 *           displayed.
 *  rstack = Array of length nrefmax containing a sequence of vertex
 *           counts (nv values) associated with refinements by 
 *           Functions refineRe and refineRt, and unrefinements by
 *           Function unrefine.
 *  sx,sy = Scale factors in the mapping from window coordinates
 *          to object coordinates -- computed in reshape.
 *  title = Optional character string written on the plot.
 *  title_font = Font number for the title string.
 *  title_pos = Raster position of (lower left corner of) title.
 *  window = Array of window identifiers.
 *  xc,yc = Object coordinates of the point that maps to the center
 *          of the viewport:  initialized to the center of the
 *          bounding box, but may be altered to pan around (when 
 *          zoomed in, for example).
 *  xmin,ymin = Coordinates of the lower left corner of the bounding
 *              box containing the vertices.
 *  xmax,ymax = Coordinates of the upper right corner of the bounding
 *              box containing the vertices.
 *  xv1,yv1 = Object coordinates of the lower left or upper right 
 *            corner of the view volume, or of an endpoint of a line
 *            segment during a drag operation:  used by Function 
 *            motion.
 *  xv2,yv2 = Object coordinates of the upper right or lower left 
 *            corner of the view volume, or of the cursor position
 *            during a drag operation:  used by Function motion.
 *  xv3,yv3 = Object coordinates of an endpoint of a line segment
 *            during a drag operation:  used by Function motion.
 */
static GLboolean antialias = GL_FALSE;
static GLboolean aspect = GL_TRUE;
static GLboolean bdry_layer = GL_FALSE;
static GLboolean color_code = GL_FALSE;
static GLfloat color_b[3] = {0.0, 0.0, 0.0};
static GLfloat color_e[3] = {1.0, 1.0, 1.0};
static GLfloat color_t[3] = {1.0, 1.0, 1.0};
static GLdouble cr = 0.0, cx, cy;
static GLboolean dlg_exit = GL_FALSE;
struct rect dlg_r0 = {240, 25, 90, 35};
struct rect dlg_r1 = {70, 25, 90, 35};
static char dlg_string[LSTRING] = "";
static GLboolean dlg_yes = GL_FALSE;
static GLdouble dx, dy;
static char *fname;
static unsigned int ind_font = 2;
static GLfloat line_width = 2.0;
static enum Mode mode = DEFAULT;
static GLint nnbv = 0;
static GLint nref = 0;
static GLint nrefmax = 10;
static GLboolean plot_box = GL_FALSE;
static GLboolean plot_e = GL_FALSE;
static GLboolean plot_nrv = GL_TRUE;
static GLboolean plot_t = GL_FALSE;
static GLboolean plot_v = GL_FALSE;
static GLint *rstack;
static GLdouble sx, sy;
static char title[LSTRING] = "";
static unsigned int title_font = 3;
static GLfloat title_pos[3] = {0., 0., 0.};
static int window[2];
static GLdouble xc, yc;
static GLdouble xmin, ymin, xmax, ymax;
static GLdouble xv1, xv2, yv1, yv2, xv3, yv3;


/*
 *  Function descriptions:
 *
 *  bbox:  Compute and store the corner coordinates of the bounding 
 *         box:  [xmin,xmax] X [ymin,ymax].
 *
 *  bcurve:  Return the CCW-ordered sequence of boundary vertex
 *           indices in a boundary curve.
 *
 *  closeR:  Close a user-selected polygon R by implicitly adding
 *           an edge from the last vertex to the first.
 *
 *  coarsenRv:  Coarsen R by removing 70% of the vertices.
 *
 *  convertRv:  Convert boundary vertices in R to removable or non-
 *              removable by deleting their indices from lnrv or 
 *              inserting their indices into lnrv if necessary.
 *
 *  deltri:  Delete a triangle.
 *
 *  delvtx:  Delete a vertex.
 *
 *  display:  Draw triangle mesh edges with annotation, title, R, etc.
 *
 *  displayPosition:  Display the coordinates of a point, such as the 
 *                    current mouse position, in the upper left corner
 *                    of the viewport.
 *
 *  displayString:  Write a string to the screen.
 *
 *  dlgDisplay:  Dialog window display callback.
 *
 *  dlgKey:  Keyboard callback associated with the dialog window.
 *
 *  dlgMouse:  Callback triggered by a mouse button press or release 
 *             with the mouse position in the dialog window.
 *
 *  dlgReshape:  Reshape function for the dialog window (window 2). 
 *
 *  dragCorner:  Mouse callback used to alter the domain by moving 
 *               a non-removable vertex.
 *
 *  drawCircle:  Draw a circle of radius cr centered at (cx,cy).
 *
 *  drawRect:  Draw the outline of an axis-aligned rectangle.
 *
 *  fillRt:  Fill R with triangles where R contains three adjacent
 *           boundary vertices forming a reentrant angle (corner) at
 *           a removable vertex.
 *
 *  getCenter:  Mouse callback used to move the center (xc,yc).
 *
 *  getCircle:  Mouse callback used to select a circle center (cx,cy)
 *              and radius cr.
 *
 *  getCorners:  Replace the set of non-removable vertices lnrv with 
 *               the set of boundary vertices that form corners (do
 *               not lie on the line defined by their predecessors
 *               and successors).
 *
 *  getLin:  Read a line from stdin.
 *
 *  getPolygon:  Mouse callback used to construct R.
 *
 *  getTitlePos:  Mouse callback used to set the title location.
 *
 *  getVertex:  Mouse callback used to select a vertex in uv edit
 *              mode.
 *
 *  getViewVolume:  Mouse callback used to change the view volume by
 *                  dragging a rectangle.
 *
 *  init1:  Initialization for window 1.
 *
 *  initR:  Initialize user-selected polygon R:  store the lists of
 *          triangles, edges, and vertices contained in R.
 *
 *  inputData:  Allocate storage and read a data set containing a tri-
 *              angulation:  nv, vtxy, nt, and first three columns of
 *              ltri.
 *
 *  insertb:  Partition a triangle into three equal-area subtriangles
 *            by inserting the barycenter as a new vertex.
 *
 *  inside:  Locate a point P relative to the polygon R, returning a
 *           nonzero value iff P is contained in R.
 *
 *  intsec:  Test for an intersection between a pair of line 
 *           segments in the plane.
 *
 *  key:  Keyboard callback for window 1 (calls menu).
 *
 *  makeMenu:  Create a menu for window 1.
 *
 *  menu:  Respond to menu selections.
 *
 *  motion:  Callback for mouse motion with a button pressed
 *
 *  mouse:  Default mouse callback:  used to alter R.
 *
 *  moveRn:  Project all non-removable vertices in R (that can be
 *           moved without destroying the mesh) onto the current
 *           circle (if any).
 *
 *  optimizeRt:  Optimize R with Delaunay edge swaps in convex    
 *               quadrilaterals consisting of pairs of adjacent 
 *               triangles in R.
 *
 *  optimizeRv:  Optimize R by applying a single iteration of a
 *               smoothing operation in which each removable vertex
 *               is moved (if possible) to a weighted sum of the 
 *               vertices in its patch (the vertex and its neighbors).
 *
 *  passiveMotion:  Callback triggered by mouse motion with no mouse
 *                  buttons pressed:  used to display the object 
 *                  coordinates of the mouse position.
 *
 *  plDistance:  Test for a point-line distance within a tolerance.
 *
 *  refineRe:  Refine R by adding a new vertex at the midpoint of each
 *             edge in R (using split and swap).
 *
 *  refineRt:  Refine R by adding triangle barycenters.
 *
 *  removeRt:  Remove the triangles in R, creating new boundary curves
 *             if necessary.
 *
 *  reshape:  Reshape callback (choose viewport and view volume)for 
 *            window 1.
 *
 *  specialKey:  Move center (xc,yc) in response to arrow keys, zoom
 *               in or out on Page Up or Page Down keys, or restore
 *               default view volume if Home key is pressed.
 *
 *  split:  Add a new vertex at the midpoint of an edge, splitting the
 *          edge into a pair of edges, and replacing each triangle 
 *          containing the edge by two triangles.
 *
 *  storeV:  Compute and store global parameter values that depend on 
 *           the input data:  dx, dy, xc, xmax, xmin, yc, ymax, and
 *           ymin, and allocate storage for R and rstack.
 *
 *  swap:  Modify a triangle mesh by swapping diagonal edges in a con-
 *         vex quadrilateral defined by a pair of adjacent triangles.
 *
 *  swapEdges:  Mouse callback used to swap edges.
 *  
 *  swptst:  Delaunay circumcircle test.
 *
 *  trmtst:  Triangulation data structure integrity test.
 *
 *  trprint:  Print triangle mesh data structure on stdout.
 *
 *  trqual1:  Compute a measure of triangle quality,
 *
 *  trqual2:  Compute a measure of triangle quality,
 *
 *  trwrite:  Output the triangle mesh to a file.
 *
 *  unrefine:  Unrefine the mesh by reversing the previous refinement
 *             by refineRe or refineRt.
 *
 *  xlate:  Translate glut window coordinates (wx,wy), with depth wz, 
 *          to object coordinates (x,y,z) by inverting the mappings 
 *          associated with the modelview, projection, and viewport 
 *          matrices.
 */


/*
 *  Function prototypes
 */
static void bbox(void);
static unsigned char bcurve(int kc, int *nbv, int indbv[]);
static void closeR(void);
static void coarsenRv(void);
static void convertRv(unsigned char removable);
static unsigned char deltri(int kt0);
static unsigned char delvtx(int kv0);
static void display(void);
static void displayPosition(GLdouble x, GLdouble y);
static int displayString(GLfloat pos[3], char *string,
                         unsigned int fontno, GLfloat scolor[3]);
static void dlgDisplay(void);
static void dlgKey(unsigned char key, int x, int y);
static void dlgMouse(int button, int state, int x, int y);
static void dlgReshape(int w, int h);
static void dragCorner(int button, int state, int x, int y);
static void drawCircle(void);
static void drawRect(GLuint x0, GLuint y0, GLuint w, GLuint h, 
                     GLfloat brcolor[3], GLfloat tlcolor[3]);
static void fillRt(void);
static void getCenter(int button, int state, int x, int y);
static void getCircle(int button, int state, int x, int y);
static void getCorners(void);
static int getLin(char *s, int len);
static void getPolygon(int button, int state, int x, int y);
static void getTitlePos(int button, int state, int x, int y);
static void getVertex(int button, int state, int x, int y);
static void getViewVolume(int button, int state, int x, int y);
static void init1(void);
static void initR(void);
static void inputData(int argc, char *argv[]);
static void insertb(int kt0);
static unsigned char inside(double xp, double yp);
static unsigned char intsec(double x1, double y1, double x2, 
                            double y2, double x3, double y3, 
                            double x4, double y4);
static void key(unsigned char key, int x, int y);
static void makeMenu(void);
static void menu(int item);
static void motion(int x, int y);
static void mouse(int button, int state, int x, int y);
static void moveRn(void);
static void optimizeRt(void);
static void optimizeRv(void);
static void passiveMotion(int x, int y);
static unsigned char plDistance(double xp, double yp, double x1, 
                                double y1, double x2, double y2, 
                                double tols);
static void refineRe(void);
static void refineRt(void);
static void removeRt(void);
static void reshape(int w, int h);
static void specialKey(int key, int x, int y);
static void split(int ke);
static void storeV(void);
static void swap(int ke, int updateR);
static void swapEdges(int button, int state, int x, int y);
static unsigned char swptst(int in1, int in2, int io1, int io2);
static void trmtst(int nc, int nb, int nv, int nn, int nt, int ne, 
                   int ltri[][9], int lvtx[], int ledg[], int lnrv[],
                   int lcrv[]);
static void trprint(int nc, int nb, int nv, int nn, int nt, int ne, 
                    double vtxy[][2], double uv[][2], int ltri[][9], 
                    int lvtx[], int ledg[], int lnrv[], int lcrv[]);
static double trqual1(double x1, double y1, double x2, double y2,
		      double x3, double y3);
static double trqual2(double x1, double y1, double x2, double y2,
		      double x3, double y3);
static void trwrite(char *fname, int nc, int nb, int nv, int nt, 
                    int ne, double vtxy[][2], double uv[][2], 
                    int ltri[][9]);
static void unrefine(void);
static void xlate(GLint wx, GLint wy, GLdouble wz, 
                  GLdouble *x, GLdouble *y, GLdouble *z);


/*
 *--------------------------------------------------------------------
 *
 *  bbox:  Compute and store the corner coordinates defining the
 *         bounding box.
 *
 *  Global variables required:  nv, vtxy, xmin, xmax, ymin, ymax
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void bbox(void)
{
   int i;
   xmin = vtxy[0][0];
   xmax = xmin;
   ymin = vtxy[0][1];
   ymax = ymin;
   for (i = 1; i < nv; i++) {
      if (vtxy[i][0] < xmin) xmin = vtxy[i][0];
      if (vtxy[i][0] > xmax) xmax = vtxy[i][0];
      if (vtxy[i][1] < ymin) ymin = vtxy[i][1];
      if (vtxy[i][1] > ymax) ymax = vtxy[i][1];
   }
   return;
}
/* End of bbox */


/*
 *--------------------------------------------------------------------
 *
 *  bcurve:  Return the CCW-ordered sequence of boundary vertex
 *           indices in a boundary curve.
 *
 *  On input:  kc is a curve index in the range 0 to nc-1.
 *
 *  On output:  The return value is 0 if kc is invalid.  Otherwise,
 *              the return value is 1, nbv is the number of vertices
 *              in the curve, and indbv contains the sequence of 
 *              indices in the first nbv positions.
 *
 *  Global variables required:  lcrv, ltri, lvtx, nc
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static unsigned char bcurve(int kc, int *nbv, int indbv[])
{
   int i, kt, kv, kv0, lbv;
/*
 *  Test for invalid input.
 */
   if (kc < 0  ||  kc >= nc) return 0;
/*
 *  Loop on boundary vertices kv starting with kv0.
 */
   lbv = 0;
   kv0 = lcrv[kc];
   kv = kv0;
   do {
      indbv[lbv] = kv;
      lbv++;
      kt = lvtx[kv];
      if (ltri[kt][0] == kv) 
         i = 1;
      else if (ltri[kt][1] == kv)
         i = 2;
      else
         i = 0;
      kv = ltri[kt][i];
   } while (kv != kv0);
   *nbv = lbv;
   return 1;
}
/* End of bcurve */


/*
 *--------------------------------------------------------------------
 *
 *  closeR:  Close the user-selected polygon R by implicitly adding
 *           an edge from the last vertex to the first.  If the added
 *           edge would result in a self-intersecting curve, the
 *           polygon is rejected by setting np = 0.
 *
 *  Global variables required:
 *
 *      np = Number of vertices in the boundary of R.
 *
 *      ner, ntr, nvr = Numbers of edges, triangles, and vertices 
 *                      in R. 
 *
 *      R = Array dimensioned npmax by 2 containing the Cartesian 
 *          coordinates of the sequence of vertices defining the 
 *          boundary of R.
 *
 *  glmesh2 function called:  intsec
 *
 *--------------------------------------------------------------------
 */
static void closeR(void)
{
   GLboolean intr = 0;
   GLint i;
/*
 *  Test for at least three boundary vertices.
 */
   if (np < 3) {
      np = 0;
      ner = 0;
      ntr = 0;
      nvr = 0;
   }
   if (np == 0) return;
/*
 *  Test the last segment for an intersection.
 */
   intr = 0;
   for (i = 1; i <= np-3; i++) {
      if (intsec(R[np-1][0],R[np-1][1],R[0][0],R[0][1],
                 R[i][0],R[i][1],R[i+1][0],R[i+1][1])) {
         intr = 1;
         break;
      }
   }
   if (intr) {
      np = 0;
      ner = 0;
      ntr = 0;
      nvr = 0;
   }
   return;
}
/* End of closeR */


/*
 *--------------------------------------------------------------------
 *
 *  coarsenRv:  Coarsen R by removing approximately .7*nvr removable
 *              vertices, chosen at random.
 *
 *  glmesh2 functions called:  initR, delvtx
 *
 *--------------------------------------------------------------------
 */
static void coarsenRv(void)
{
   unsigned char removable;
   int i, k, kv, n, nd, nrv;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
/*
 *  Count the number nrv of removable vertices in R, and store
 *  their indices in lste.
 */
   nrv = 0;
   for (k = 0; k < nvr; k++) {
      kv = listvr[k];
      removable = 1;
      for (i = 0; i < nn; i++) {
         if (lnrv[i] == kv) {
            removable = 0;
            break;
         }
      } 
      if (!removable) continue;
      lste[nrv] = kv;
      nrv++;
   }
/*
 *  Delete n = max(.7*nvr,nrv) removable vertices.
 */
   n = (int) (0.7*(double) nvr + 0.5);
   if (n < nrv) n = nrv;
   nd = 0;
   while (nd < n  &&  nrv > 0) {
      k = rand()%nrv;
      kv = lste[k];
      if (delvtx(kv)) nd++;
      nrv--;
      lste[kv] = lste[nrv];
   }
/*
 *  The refinement stack is emptied if a vertex was deleted.
 */
   if (nd) nref = 0;
   return;
}
/* End of coarsenRv */


/*
 *--------------------------------------------------------------------
 *
 *  convertRv:  Convert boundary vertices in R to removable or non-
 *              removable by deleting their indices from lnrv or 
 *              inserting their indices into lnrv if necessary.
 *
 *  On input:  removable is a flag specifying that the boundary
 *             vertices in R are to be converted to removable (TRUE)
 *             or non-removable (FALSE).
 *
 *  glmesh2 function called:  initR
 *
 *--------------------------------------------------------------------
 */
static void convertRv(unsigned char removable)
{
   unsigned char found;
   int i, j, k, kt, kv;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
   if (removable) {
/*
 *  For each boundary vertex kv in R, search lnrv for kv, and 
 *  remove kv if found.
 */
      for (k = 0; k < nvr; k++) {
         kv = listvr[k];
         kt = lvtx[kv];
         if (ltri[kt][0] == kv)
            j = 2;
         else if (ltri[kt][1] == kv)
            j = 0;
         else
            j = 1;
         if (ltri[kt][j+3] >= 0) continue; 
         for (i = 0; i < nn; i++) {
            if (lnrv[i] == kv) { 
               nn--;
               lnrv[i] = lnrv[nn];
               break;
            }
         }
      }
   } else {
/*
 *  For each boundary vertex kv in R, search lnrv for kv, and add kv
 *  if not found.
 */
      for (k = 0; k < nvr; k++) {
         kv = listvr[k];
         kt = lvtx[kv];
         if (ltri[kt][0] == kv)
            j = 2;
         else if (ltri[kt][1] == kv)
            j = 0;
         else
            j = 1;
         if (ltri[kt][j+3] >= 0) continue; 
         found = 0;
         for (i = 0; i < nn; i++) {
            if (lnrv[i] == kv) { 
               found = 1;
               break;
            }
         }
         if (!found) {
            lnrv[nn] = kv;
            nn++;
         }
      }
   }
   return;
}
/* End of convertRv */


/*
 *--------------------------------------------------------------------
 *
 *  deltri:  Delete a triangle.
 *
 *  This function removes a triangle from the triangle mesh unless its
 *  removal would delete a non-removable vertex or destroy the 
 *  mesh by leaving the union of remaining triangles not strongly 
 *  connected.
 *
 *  On input:  kt0 is the index of the triangle to be removed.
 *
 *  On output:  the triangle is removed unless its index is not valid
 *              (in the range 0 to nt-1), or it is not in R, or its 
 *              removal would render the union of triangles weakly 
 *              connected.  The return value is the number of 
 *              triangles removed.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static unsigned char deltri(int kt0)
{
   int i, j, k, ke, kt, kt1, ktn, kv, kvs, nbe;

/*
 *  Test for invalid index.
 */
   if (kt0 < 0  ||  kt0 >= nt) return 0;
/*
 *  kt0 cannot be removed if it has a boundary vertex kv which is
 *  not an endpoint of a boundary edge contained in kt0.  Also,
 *  test for a vertex kv of kt0 not in R.
 */
   for (i = 0; i < 3; i++) {
      kv = ltri[kt0][i];
      if (!listf[kv]) return 0;
      kt = lvtx[kv];
      if (ltri[kt][0] == kv)
         j = 2;
      else if (ltri[kt][1] == kv)
         j = 0;
      else 
         j = 1;
      if (ltri[kt][j+3] < 0  &&  kt != kt0  &&  
          ltri[kt0][(i+1)%3+3] >= 0) return 0;
   }
/*
 *  Test for a non-removable vertex shared by two boundary edges
 *  in kt0.
 */
   for (i = 0; i < 3; i++) {
      kt = ltri[kt0][(i+1)%3+3];
      ktn = ltri[kt0][(i+2)%3+3];
      if (kt < 0  &&  ktn < 0) {
         kv = ltri[kt0][i];
         for (k = 0; k < nn; k++) {
            if (lnrv[k] == kv) return 0;
         }
      }
   }
/*
 *  Save a vertex index for updating lcrv if necessary.
 */
   kvs = kv;
/*
 *  Loop on vertices kv and opposite edges ke of triangle kt0,
 *  accumulating nbe (number of boundary edges in kt0).
 */
   nbe = 0;
   for (i = 0; i < 3; i++) {
      kv = ltri[kt0][i];
/*
 *  Alter lvtx[kv] if necessary, to reflect the new boundary curve.
 */
      kt = ltri[kt0][(i+1)%3+3];
      if (kt >= 0) {
         lvtx[kv] = kt;
      } else {
         if (ltri[kt0][(i+2)%3+3] < 0) {
/*
 *  Remove vertex kv, and relabel vertex nv-1 to kv.
 */            
            nv--;
            vtxy[kv][0] = vtxy[nv][0];
            vtxy[kv][1] = vtxy[nv][1];
            uv[kv][0] = uv[nv][0];
            uv[kv][1] = uv[nv][1];
            kt = lvtx[nv];
            lvtx[kv] = kt;
            kt1 = kt;
            do {
               if (ltri[kt][0] == nv)
                  j = 0;
               else if (ltri[kt][1] == nv)
                  j = 1;
               else
                  j = 2;
               ltri[kt][j] = kv;
               j = (j+1)%3;
               kt = ltri[kt][j+3];
            } while (kt >= 0  &&  kt != kt1);
/*
 *  Correct lcrv if necessary.
 */
            for (j = 0; j < nc; j++) {
               if (lcrv[j] == kv) lcrv[j] = ltri[kt0][(i+1)%3];
               if (lcrv[j] == nv) lcrv[j] = kv;
            }
/*
 *  Correct lnrv if necessary.
 */
            for (j = 0; j < nn; j++) {
               if (lnrv[j] == nv) {
                  lnrv[j] = kv;
                  break;
               }
            }
/*
 *  Remove kv from listvr.
 */
            for (k = 0; k < nvr; k++) {
               if (listvr[k] == kv) {
                  nvr--;
                  listvr[k] = listvr[nvr];
                  break;
               }
            }
            for (k = 0; k < nvr; k++) {
               if (listvr[k] == nv) {
                  listvr[k] = kv;
                  break;
               }
            }
/*
 *  Adjust listf.
 */
            listf[kv] = listf[nv];
         }
      }
      ktn = ltri[kt0][i+3];
      ke = ltri[kt0][i+6];
      if (ktn >= 0) {
/*
 *  Remove kt0 as a neighbor of ktn, and alter ledg[ke] if necessary.
 */
         if (ltri[ktn][3] == kt0)
            j = 0;
         else if (ltri[ktn][4] == kt0)
            j = 1;
         else
            j = 2;
         ltri[ktn][j+3] = -1;
         ledg[ke] = ktn;
      } else {
/*
 *  Remove boundary edge ke, and relabel edge ne-1 to ke.
 */
         nbe++;
         ne--;
         kt = ledg[ne];
         ledg[ke] = kt;
         if (ltri[kt][6] == ne)
            j = 6;
         else if (ltri[kt][7] == ne)
            j = 7;
         else
            j = 8;
         ltri[kt][j] = ke;
         kt = ltri[kt][j-3];
         if (kt >= 0) {
            if (ltri[kt][6] == ne)
               j = 6;
            else if (ltri[kt][7] == ne)
               j = 7;
            else
               j = 8;
            ltri[kt][j] = ke;
         }
/*
 *  Remove ke from lister.
 */
         for (k = 0; k < ner; k++) {
            if (lister[k] == ke) {
               ner--;
               lister[k] = lister[ner];
               break;
            }
         }
         for (k = 0; k < ner; k++) {
            if (lister[k] == ne) {
               lister[k] = ke;
               break;
            }
         }
      }
   }
/*
 *  Remove triangle kt0, and relabel triangle nt-1 to kt0.
 */
   nt--;
   if (kt0 != nt) {
      for (j = 0; j < 9; j++) {
         ltri[kt0][j] = ltri[nt][j];
      }
      for (i = 0; i < 3; i++) {
         kv = ltri[nt][i];
         if (lvtx[kv] == nt) lvtx[kv] = kt0;
         kt = ltri[nt][i+3];
         if (kt >= 0) {
            if (ltri[kt][3] == nt)
               j = 3;
            else if (ltri[kt][4] == nt)
               j = 4;
            else 
               j = 5;
            ltri[kt][j] = kt0;
         }
         ke = ltri[nt][i+6];
         if (ledg[ke] == nt) {
            if (kt > kt0) { 
               ledg[ke] = kt;
            } else {
               ledg[ke] = kt0;
            }
         }
      }
   }
/*
 *  Remove kt0 from listtr.
 */
   for (k = 0; k < ntr; k++) {
      if (listtr[k] == kt0) {
         ntr--;
         listtr[k] = listtr[ntr];
         break;
      }
   }
   for (k = 0; k < ntr; k++) {
      if (listtr[k] == nt) {
         listtr[k] = kt0;
         break;
      }
   }
/*
 *  Adjust counts nb and nc, and update lcrv if a new curve was added.
 */
   if (nbe == 0) {
      lcrv[nc] = kvs;
      nc++;
      nb += 3;
   } else if (nbe == 1) {
      nb++;
   } else if (nbe == 2) {
      nb--;
   } else {
      nc--;
      nb -= 3;
   }
   return 1;
}
/* End of deltri */


/*
 *--------------------------------------------------------------------
 *
 *  delvtx:  Delete a vertex.
 *
 *  All possible edge swaps (Function swap) are applied to the edges
 *  incident on vertex kv0 so that kv0, if an interior vertex,  has 
 *  three neighbors, and then kv0 and the edges containing kv0 are
 *  removed, replacing three triangles by two triangles in the case of
 *  an interior vertex.  The vertex is thus removed without leaving a
 *  hole.
 *
 *  On input:  kv0 is the index of the vertex to be removed.
 *
 *  On output:  The return value is the number of vertices removed
 *              (0 or 1).  If index kv0 is invalid, or kv0 is an
 *              element of lnrv, no variables are altered.  It is
 *              possible that kv0 indexes an interior vertex but the
 *              number of neighbors of kv0 could not be reduced to
 *              three due to collinear vertices and floating point
 *              inaccuracy.  Also, it is possible that a boundary 
 *              vertex cannot be removed without destroying the 
 *              triangle mesh.  The mesh may be altered by swaps but
 *              kv0 is not removed in these cases.
 *
 *  Note that vertex kv0 may be deleted even if it is not contained 
 *  in R.
 *
 *  glmesh2 function called:  swap
 *
 *--------------------------------------------------------------------
 */
static unsigned char delvtx(int kv0)
{
   unsigned char bndry, swp;
   int i, i0, iter, j, k, ke, ke0, ken, kt, kt0, kt1, ktn, kv, kv1, 
       kv2, kv3, kvs = 0, nnb;
/*
 *  Test for invalid input.
 */
   if (kv0 < 0  ||  kv0 >= nv  ||  nv < 4) return 0;
/*
 *  Test for kv0 in lnrv.
 */
   for (i = 0; i < nn; i++) {
      if (lnrv[i] == kv0) return 0;
   }
/*
 *  Count the number of neighbors of kv0.
 *
 *  kt1 = index of first triangle containing kv0.
 *  nnb = number of neighbors of kv0.
 */
   kt1 = lvtx[kv0];
   kt = kt1;
   nnb = 0;
   do {
      if (ltri[kt][0] == kv0)
         i = 1;
      else if (ltri[kt][1] == kv0)
         i = 2;
      else
         i = 0;
      kt = ltri[kt][i+3];
      nnb++;
   } while (kt >= 0  &&  kt != kt1);
   bndry = (kt < 0);
   if (bndry) nnb++;
/*
 *  Loop on edges ke incident on kv0 and shared by triangles
 *  kt and ktn, swapping out ke if possible.
 *
 *  bndry = 1 iff kv0 is a boundary vertex.
 *  kt1 = index of first triangle containing kv0.
 *  nnb = number of neighbors of kv0.
 *  swp = flag with value 1 initially and following a swap.
 */
   kt = kt1;
   swp = 1;
   while ((swp  ||  kt != kt1)  &&  (bndry  ||  nnb > 3)) {
      if (ltri[kt][0] == kv0)
         i = 1;
      else if (ltri[kt][1] == kv0)
         i = 2;
      else
         i = 0;
      ktn = ltri[kt][i+3];
      if (ktn < 0) break;
      ke = ltri[kt][i+6];
      kv2 = ltri[kt][i];
      kv1 = ltri[kt][(i+1)%3];
      if (ltri[ktn][0] == kv0)
         i = 2;
      else if (ltri[ktn][1] == kv0)
         i = 0;
      else
         i = 1;
      kv3 = ltri[ktn][i];
/*
 *  Edge ke can be swapped iff kv1 strictly Left kv3 -> kv2 and 
 *  kv0 Left kv2 -> kv3.
 */
      swp = ((vtxy[kv2][0]-vtxy[kv3][0])*(vtxy[kv1][1]-vtxy[kv3][1]) >
             (vtxy[kv2][1]-vtxy[kv3][1])*(vtxy[kv1][0]-vtxy[kv3][0])       
          && (vtxy[kv3][0]-vtxy[kv2][0])*(vtxy[kv0][1]-vtxy[kv2][1]) >=
             (vtxy[kv3][1]-vtxy[kv2][1])*(vtxy[kv0][0]-vtxy[kv2][0]));
      if (swp) {
         swap(ke, 1); 
         nnb--;
         if (kt == kt1) kt1 = ktn;
      }
      kt = ktn;
   }
/*
 *  Test for failure with kv0 interior and having more than three 
 *  neighbors, or kv0 a boundary vertex with a boundary vertex as
 *  a neighbor other than the first or last.
 */
   if (!bndry) {
      if (nnb > 3) return 0;
   } else {
      kt = lvtx[kv0];
      while (kt >= 0) {
         if (ltri[kt][0] == kv0)
            i = 1;
         else if (ltri[kt][1] == kv0)
            i = 2;
         else
            i = 0;
         ktn = ltri[kt][i+3];
         if (ktn < 0) break;
         i = (i+1)%3;
         kv1 = ltri[kt][i];
         kt1 = lvtx[kv1];
         if (ltri[kt1][0] == kv1)
            j = 2;
         else if (ltri[kt1][1] == kv1)
            j = 0;
         else
            j = 1;
         if (ltri[kt1][j+3] < 0) return 0;
         kt = ktn;
      }
   }
/*
 *  Set kt0 and ke0 to the indices of the first triangle and edge,
 *  respectively, to be deleted, and set i0 to the ltri column
 *  index of kv0 in kt1.
 */
   kt1 = lvtx[kv0];
   if (ltri[kt1][0] == kv0)
      i0 = 0;
   else if (ltri[kt1][1] == kv0)
      i0 = 1;
   else
      i0 = 2;
   i = i0;
   if (!bndry) {
      i = (i+1)%3;
      ke0 = ltri[kt1][i+6];
      kt0 = ltri[kt1][i+3];
      j = (i+1)%3;
      kv1 = ltri[kt1][j];
      if (lvtx[kv1] == kt0) lvtx[kv1] = kt1;
   } else {
/* 
 *  Save a neighbor kvs of kv0 for updating lcrv if necessary.
 */
      kvs = ltri[kt1][(i0+1)%3];
      i = (i+2)%3;
      ke0 = ltri[kt1][i+6];
      kt0 = kt1;
      kt1 = -1;
   }
/*
 *  Remove edge ke0, and relabel edge ne-1 to ke0.
 */
   ne--;
   if (ke0 != ne) {
      kt = ledg[ne];
      ledg[ke0] = kt;
      if (ltri[kt][6] == ne)
         j = 6;
      else if (ltri[kt][7] == ne)
         j = 7;
      else
         j = 8;
      ltri[kt][j] = ke0;
      kt = ltri[kt][j-3];
      if (kt >= 0) {
         if (ltri[kt][6] == ne)
            j = 6;
         else if (ltri[kt][7] == ne)
            j = 7;
         else
            j = 8;
         ltri[kt][j] = ke0;
      }
   }
/*
 *  Remove ke0 from lister.
 */
   for (k = 0; k < ner; k++) {
      if (lister[k] == ke0) {
         ner--;
         lister[k] = lister[ner];
         break;
      }
   }
   for (k = 0; k < ner; k++) {
      if (lister[k] == ne) {
         lister[k] = ke0;
         break;
      }
   }
/*
 *  Loop on triangles kt0 with triangle ktn and edge ken opposite
 *  vertex kv0.  Note that kt1 = -1 iff bndry = 1.
 */
   for (iter = 1; iter < nnb; iter++) {
      if (ltri[kt0][0] == kv0)
         i = 0;
      else if (ltri[kt0][1] == kv0)
         i = 1;
      else
         i = 2;
      ktn = ltri[kt0][i+3];
      ken = ltri[kt0][i+6];
      if (ktn >= 0) {
         if (ltri[ktn][3] == kt0)
            ltri[ktn][3] = kt1;
         else if (ltri[ktn][4] == kt0)
            ltri[ktn][4] = kt1;
         else
            ltri[ktn][5] = kt1;
      }
      if (!bndry) {
         j = (i+1)%3;
         kv1 = ltri[kt0][j];
         j = (i+2)%3;
         kv2 = ltri[kt0][j];
         if (lvtx[kv1] == kt0) lvtx[kv1] = kt1;
         if (lvtx[kv2] == kt0) lvtx[kv2] = kt1;
         ledg[ken] = (kt1 > ktn) ? kt1 : ktn;
         ltri[kt1][(i0+iter)%3+3] = ktn;
         ltri[kt1][(i0+iter)%3+6] = ken;
      } else {
         j = (i+2)%3;
         kv2 = ltri[kt0][j];
         lvtx[kv2] = ktn;
         ledg[ken] = ktn;
      }
/*
 *  Set ke0 to the edge of triangle kt0 that contains kv0 and was not
 *  removed at the previous iteration (or before the loop), and set
 *  ktn to the next triangle to be processed.
 */
      i = (i+1)%3;
      ke0 = ltri[kt0][i+6];
      ktn = ltri[kt0][i+3];
/*
 *  Remove edge ke0, and relabel edge ne-1 to ke0.
 */
      ne--;
      if (ke0 != ne) {
         kt = ledg[ne];
         ledg[ke0] = kt;
         if (ltri[kt][6] == ne)
            j = 6;
         else if (ltri[kt][7] == ne)
            j = 7;
         else
            j = 8;
         ltri[kt][j] = ke0;
         kt = ltri[kt][j-3];
         if (kt >= 0) {
            if (ltri[kt][6] == ne)
               j = 6;
            else if (ltri[kt][7] == ne)
               j = 7;
            else
               j = 8;
            ltri[kt][j] = ke0;
         }
      }
/*
 *  Remove ke0 from lister.
 */
      for (k = 0; k < ner; k++) {
         if (lister[k] == ke0) {
            ner--;
            lister[k] = lister[ner];
            break;
         }
      }
      for (k = 0; k < ner; k++) {
         if (lister[k] == ne) {
            lister[k] = ke0;
            break;
         }
      }
/*
 *  Store the new first vertex of kt1 if kv0 is interior.
 */        
      if (!bndry  &&  iter == 1) ltri[kt1][i0] = ltri[kt0][(i+1)%3];
/*
 *  Remove triangle kt0, and relabel triangle nt-1 to kt0.
 */
      nt--;
      if (kt0 != nt) {
         if (kt1 == nt) kt1 = kt0;
         if (ktn == nt) ktn = kt0;
         for (j = 0; j < 9; j++) {
            ltri[kt0][j] = ltri[nt][j];
         }
         for (i = 0; i < 3; i++) {
            kv = ltri[nt][i];
            if (lvtx[kv] == nt) lvtx[kv] = kt0;
            kt = ltri[nt][i+3];
            if (kt >= 0) {
               if (ltri[kt][3] == nt)
                  ltri[kt][3] = kt0;
               else if (ltri[kt][4] == nt)
                  ltri[kt][4] = kt0;
               else if (ltri[kt][5] == nt)
                  ltri[kt][5] = kt0;
            }
            ke = ltri[nt][i+6];
            if (ledg[ke] == nt) {
               ledg[ke] = (kt > kt0) ? kt : kt0;
            }
         }
      }
/*
 *  Remove kt0 from listtr.
 */
      for (k = 0; k < ntr; k++) {
         if (listtr[k] == kt0) {
            ntr--;
            listtr[k] = listtr[ntr];
            break;
         }
      }
      for (k = 0; k < ntr; k++) {
         if (listtr[k] == nt) {
            listtr[k] = kt0;
            break;
         }
      }
/*
 *  Bottom of for loop.
 */
      kt0 = ktn;
   }
/*
 *  Remove vertex kv0 from vtxy, uv, and lvtx, and relabel vertex 
 *  nv-1 to kv0.
 */            
   nv--;
   if (bndry) nb += nnb-3;
   for (j = 0; j < nc; j++) {
      if (lcrv[j] == kv0) lcrv[j] = kvs;
      if (lcrv[j] == nv) lcrv[j] = kv0;
   }
   if (kv0 != nv) {
      vtxy[kv0][0] = vtxy[nv][0];
      vtxy[kv0][1] = vtxy[nv][1];
      uv[kv0][0] = uv[nv][0];
      uv[kv0][1] = uv[nv][1];
      kt = lvtx[nv];
      lvtx[kv0] = kt;
      kt1 = kt;
      do {
         if (ltri[kt][0] == nv)
            j = 0;
         else if (ltri[kt][1] == nv)
            j = 1;
         else
            j = 2;
         ltri[kt][j] = kv0;
         j = (j+1)%3;
         kt = ltri[kt][j+3];
      } while (kt >= 0  &&  kt != kt1);
   }
/*
 *  Remove kv0 from listvr.
 */
   for (k = 0; k < nvr; k++) {
      if (listvr[k] == kv0) {
         nvr--;
         listvr[k] = listvr[nvr];
         break;
      }
   }
   for (k = 0; k < nvr; k++) {
      if (listvr[k] == nv) {
         listvr[k] = kv0;
         break;
      }
   }
/*
 *  Relabel nv in lnrv.
 */
   for (i = 0; i < nn; i++) {
      if (lnrv[i] == nv) { 
         lnrv[i] = kv0;
         break;
      }
   }
/*
 *  Adjust listf.
 */
   listf[kv0] = listf[nv];
   return 1;
}
/* End of delvtx */


/*
 *--------------------------------------------------------------------
 *
 *  dlgDisplay:  Dialog window display callback.
 *
 *  Global variables required:
 *
 *      dlg_r0 = Rectangle representing 'No' button.
 *
 *      dlg_r1 = Rectangle representing 'Yes' button.
 *
 *      dlg_string = Text string displayed in dialog box window.
 *
 *  glmesh2 functions called:  displayString, drawRect
 *
 *--------------------------------------------------------------------
 */
static void dlgDisplay(void)
{
   GLfloat black[3] = {0.0, 0.0, 0.0};   /* Border colors */
   GLfloat white[3] = {1.0, 1.0, 1.0};
   int fontno = 3;    /* Font number for strings (24-point
                           proportional spaced Times-Roman) */
   GLfloat pos[3];    /* Text string position */
   char text[7];      /* Text string */
   int x0, y0;        /* Coordinates of the lower left corner of a 
                         rectangle */
/*
 *  Clear the color buffer to the background color.
 */
   glClear(GL_COLOR_BUFFER_BIT);
/*
 *  Display dlg_string.
 */
   pos[0] = 75.0f;
   pos[1] = 160.0f;
   pos[2] = 0.0f;
   displayString(pos, dlg_string, fontno, white);
/*
 *  Draw buttons at the bottom of the window.
 */
   x0 = dlg_r1.x0;
   y0 = dlg_r1.y0;
   drawRect(x0, y0, dlg_r1.w, dlg_r1.h, black, white);
   pos[0] = (GLfloat) (x0+7);
   pos[1] = (GLfloat) (y0+7);
   sprintf(text, "  YES ");
   displayString(pos, text, fontno, white);

   x0 = dlg_r0.x0;
   y0 = dlg_r0.y0;
   drawRect(x0, y0, dlg_r0.w, dlg_r0.h, black, white);
   pos[0] = (GLfloat) (x0+7);
   pos[1] = (GLfloat) (y0+7);
   sprintf(text, "  NO  ");
   displayString(pos, text, fontno, white);

   glFlush();
   return;
}
/* End of dlgDisplay */


/*
 *--------------------------------------------------------------------
 *
 *  dlgKey:  Keyboard callback associated with the dialog window:
 *              Enter ==> Yes
 *              Escape ==> No
 *
 *  The mouse coordinates (x,y) are not used.
 *
 *  Global variables required:
 *
 *      dlg_r1 = Rectangle representing 'Yes' button.
 *
 *      dlg_yes = Flag with value defined by the 'Yes' button:  set
 *                to FALSE by the Escape key.
 *
 *  glmesh2 function called:  dlgMouse
 *
 *--------------------------------------------------------------------
 */
static void dlgKey(unsigned char key, int x, int y)
{
   if (key == 10  ||  key == 13) {        /* Enter key pressed */
      dlgMouse(GLUT_LEFT_BUTTON, GLUT_UP,
               dlg_r1.x0+1, 198-dlg_r1.y0);
   } else if (key == 27) {                /* Escape key pressed */
      dlg_yes = GL_FALSE;
      glutHideWindow();
      glutSetWindow(window[0]);
      glutShowWindow();          /* Restore visibility of window 1 */
   }
   return;
}
/* End of dlgKey */


/*
 *--------------------------------------------------------------------
 *
 *  dlgMouse:  Callback triggered by a mouse button press or release 
 *             with the mouse position in the dialog window.
 *
 *  If the left mouse button is pressed with the mouse position 
 *  inside the 'Yes' or 'No' rectangle, the appearance of the 
 *  rectangle is altered.  If the left button is released inside one 
 *  of these rectangles, dlg_yes is set to TRUE or FALSE, and the 
 *  window state is changed to hidden.
 *
 *  Global variables required:
 *
 *      dlg_exit = Flag with value TRUE iff the program is to be
 *                 terminated when the user selects the 'Yes' button.
 *
 *      dlg_r0 = Rectangle representing 'No' button.
 *
 *      dlg_r1 = Rectangle representing 'Yes' button.
 *
 *      dlg_yes = Flag with value TRUE iff the user selects the 'Yes'
 *                button in the dialog box.
 *
 *      fname = File name for output.
 *
 *      Triangulation variables:  ltri, nb, nc, ne, nt, nv, vtxy, uv
 *
 *  glmesh2 functions called:  drawRect, trwrite
 *
 *--------------------------------------------------------------------
 */
static void dlgMouse(int button, int state, int x, int y)
{
   GLfloat black[3] = {0.0, 0.0, 0.0};   /* Border colors */
   GLfloat white[3] = {1.0, 1.0, 1.0};
   int x0, y0;   /* Lower left corner of a rectangle */
   int x1, y1;   /* Upper right corner of a rectangle */

   y = 199 - y;  /* Map glut y-coordinate to viewport */

/*
 *  Left button press or release.
 */
   x0 = dlg_r1.x0;
   y0 = dlg_r1.y0;
   x1 = x0 + dlg_r1.w - 1;
   y1 = y0 + dlg_r1.h - 1;
   if (x > x0  &&  x < x1  &&  y > y0  &&  y < y1) {
/*
 *  'Yes' button.
 */
      if (state == GLUT_DOWN) {
         drawRect(x0, y0, dlg_r1.w, dlg_r1.h, white, black);
         glFlush();
      } else {                         /* 'Yes' button released */
         if (dlg_exit) {
/*
 *  Write the triangle mesh to a file, and terminate execution.
 */
            trwrite(fname, nc, nb, nv, nt, ne, vtxy, uv, ltri);
            exit(0);
         }
         dlg_yes = GL_TRUE;
         glutHideWindow();
         glutSetWindow(window[0]);
         glutShowWindow();
      }
      return;
   }

   x0 = dlg_r0.x0;
   y0 = dlg_r0.y0;
   x1 = x0 + dlg_r0.w - 1;
   y1 = y0 + dlg_r0.h - 1;
   if (x > x0  &&  x < x1  &&  y > y0  &&  y < y1) {
/*
 *  'No' button.
 */
      if (state == GLUT_DOWN) {
         drawRect(x0, y0, dlg_r0.w, dlg_r0.h, white, black);
         glFlush();
      } else {                /* 'No' button released */
         dlg_yes = GL_FALSE;
         glutHideWindow();
         glutSetWindow(window[0]);
         glutShowWindow(); 
      }
      return;
   }
   return;
}
/* End of dlgMouse */


/*
 *--------------------------------------------------------------------
 *
 *  dlgReshape:  Set the viewport, projection matrix, and modelview 
 *               matrix for the dialog window (window 2):  
 *               two-dimensional rendering with a 400 by 200 viewport.
 *
 *  This function is called before the first call to dlgDisplay.
 *
 *  If the window is too small, a request for a larger window is 
 *  executed.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void dlgReshape(int w, int h)
{
   if (w < 400  ||  h < 200) {       /* Test for sufficient */
      glutReshapeWindow(400, 200);   /*   window extent.    */
      return;
   }
/*
 *  Set the viewport to the entire window (the default).
 */
   glViewport(0, 0, (GLint) w, (GLint) h);
/*
 *  Set up an orthogonal projection for data space equal to
 *  viewport (identity).
 */
   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
   gluOrtho2D(0.0, (GLdouble) w, 0.0, (GLdouble) h);
/*
 *  Set the modelview matrix to the identity.
 */
   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();

   return;
}
/* End of dlgReshape */


/*
 *--------------------------------------------------------------------
 *
 *  display:  Clear the color buffer and draw the triangle mesh edges
 *            with optional index annotation of vertices, edges, and
 *            triangles.  The title, and selected polygon R are also 
 *            drawn, and a bounding box may be rendered.
 *
 *  glmesh2 functions called:  displayString, drawCircle, trqual1
 *
 *--------------------------------------------------------------------
 */
static void display(void)
{
   unsigned int fontno;        /* Font for labels */
   GLfloat pos[3];             /* Label position */
   char text[30];              /* Label text */

   GLfloat red[3] = {1.0, 0.0, 0.0};
   GLfloat green[3] = {0.0, 1.0, 0.0};
   GLfloat blue[3] = {0.0, 0.0, 1.0};
   GLfloat grey[3] = {0.7, 0.7, 0.7};
   double q;
   int i;
   GLdouble ch;
   GLint ke, kt, kv, kv1, kv2, kv3;
/*
 *  Clear the color buffer, and set line and point widths.
 */
   glClear(GL_COLOR_BUFFER_BIT);
   glLineWidth(line_width);
   glPointSize(line_width+1);
   if (color_code) {
/*
 *  Render constant-color triangles using color-coding based on
 *  quality.
 */
      glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
      for (kt = 0; kt < nt; kt++) {
	 kv1 = ltri[kt][0];
	 kv2 = ltri[kt][1];
	 kv3 = ltri[kt][2];
         q = trqual1(vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0],
                     vtxy[kv2][1], vtxy[kv3][0], vtxy[kv3][1]);
         if (q < 0.600)
            glColor3fv(red);
         else if (q < 0.866) 
            glColor3fv(blue);
         else
            glColor3fv(green);
	 glBegin(GL_TRIANGLES);
	 glVertex2d(vtxy[kv1][0], vtxy[kv1][1]);
	 glVertex2d(vtxy[kv2][0], vtxy[kv2][1]);
	 glVertex2d(vtxy[kv3][0], vtxy[kv3][1]);
         glEnd();
      }
   }
/*
 *  Draw the triangle mesh edges.
 */
   glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
   glColor3fv(color_e);   
   for (ke = 0; ke < ne; ke++) {
      kt = ledg[ke];
      if (ltri[kt][6] == ke)
         i = 1;
      else if (ltri[kt][7] == ke)
         i = 2;
      else
         i = 0;
      kv1 = ltri[kt][i];
      kv2 = ltri[kt][(i+1)%3];
      glBegin(GL_LINES);
         glVertex2d(vtxy[kv1][0], vtxy[kv1][1]);
         glVertex2d(vtxy[kv2][0], vtxy[kv2][1]);
      glEnd();
   }
/*
 *  Display the selected polygon R (both line segments and vertices).
 */
   glColor3fv(red);
   for (i = 0; i < np; i++) {
      glBegin(GL_LINES);
         glVertex2d(R[i][0], R[i][1]);
         glVertex2d(R[(i+1)%np][0], R[(i+1)%np][1]);
      glEnd();
      glBegin(GL_POINTS);
         glVertex2d(R[i][0], R[i][1]);
      glEnd();
   }
   if (cr > 0.0) {
/*
 *  Draw a circle of radius cr centered at (cx,cy).
 */
      glColor3fv(grey);
      drawCircle();
   }
   if (plot_v) {
/*
 *  Display vertex indices.
 */
      fontno = ind_font;
      pos[2] = 0.0;
      for (kv = 0; kv < nv; kv++) {
         pos[0] = (GLfloat) (vtxy[kv][0] + .005*dx);
         pos[1] = (GLfloat) (vtxy[kv][1] + .005*dy);
         sprintf(text, "%u", kv);
         displayString(pos, text, fontno, color_t);
      }
   }
   if (plot_nrv) {
/*
 *  Display non-removable vertex indices.
 */
      fontno = ind_font;
      pos[2] = 0.0;
      for (i = 0; i < nn; i++) {
         kv = lnrv[i];
         pos[0] = (GLfloat) (vtxy[kv][0] + .005*dx);
         pos[1] = (GLfloat) (vtxy[kv][1] + .005*dy);
         sprintf(text, "%u", kv);
         displayString(pos, text, fontno, green);
      }
   }
   if (plot_e) {
/*
 *  Display edge indices.
 */
      fontno = ind_font;
      pos[2] = 0.0;
      for (ke = 0; ke < ne; ke++) {
         kt = ledg[ke];
         if (ltri[kt][6] == ke)
            i = 1;
         else if (ltri[kt][7] == ke)
            i = 2;
         else
            i = 0;
         kv1 = ltri[kt][i];
         kv2 = ltri[kt][(i+1)%3];
         pos[0] = (GLfloat) ((vtxy[kv1][0]+vtxy[kv2][0])/2.0
                             + .005*dx);
         pos[1] = (GLfloat) ((vtxy[kv1][1]+vtxy[kv2][1])/2.0
                             + .005*dy);
         sprintf(text, "%u", ke);
         displayString(pos, text, fontno, color_t);
      }
   }
   if (plot_t) {
/*
 *  Display triangle indices.
 */
      fontno = ind_font;
      pos[2] = 0.0;
      for (kt = 0; kt < nt; kt++) {
         kv1 = ltri[kt][0];
         kv2 = ltri[kt][1];
         kv3 = ltri[kt][2];
         pos[0] = (GLfloat) ((vtxy[kv1][0]+vtxy[kv2][0]+vtxy[kv3][0])
                             /3.0);
         pos[1] = (GLfloat) ((vtxy[kv1][1]+vtxy[kv2][1]+vtxy[kv3][1])
                             /3.0);
         sprintf(text, "%u", kt);
         displayString(pos, text, fontno, color_t);
      }
   }
   if (plot_box) {
/*
 *  Display the bounding box with the lower left and upper right 
 *  corners labeled.
 */
      glColor3fv(grey);
      fontno = 5;
      glBegin(GL_LINE_LOOP);
        glVertex2d(xmin, ymin);
        glVertex2d(xmax, ymin);
        glVertex2d(xmax, ymax);
        glVertex2d(xmin, ymax);
      glEnd();
/*
 *  The character height ch is approximately 11 pixels (converted
 *  to object coordinates).
 */
      ch = (sx <= sy) ? sx : sy;   /* ch = min(sx, sy) */
      ch = 11.0*ch;
      pos[0] = (GLfloat) (xmin);
      pos[1] = (GLfloat) (ymin - 1.5*ch);
      pos[2] = 0.0;
      sprintf(text, "(%3.2f, %3.2f)", xmin, ymin);
      displayString(pos, text, fontno, color_t);
      pos[0] = (GLfloat) (xmax - 6.0*ch);
      pos[1] = (GLfloat) (ymax + .5*ch);
      sprintf(text, "(%3.2f, %3.2f)", xmax, ymax);
      displayString(pos, text, fontno, color_t);
   }

   if (*title != '\0') {
/*
 *  Display the title string
 */
      displayString(title_pos, title, title_font, color_t);
   }

   glColor3fv(color_e); 
   glutSwapBuffers();
   return;
}
/* End of display */


/*
 *--------------------------------------------------------------------
 *
 *  displayPosition:  Display the coordinates of a point, such as the 
 *                    current mouse position, in the upper left corner
 *                    of the viewport.
 *
 *  On input:  x and y are the coordinates to be displayed.
 *
 *  Global variables required:  color_b, color_t, dx, dy, sx, sy, 
 *                              xc, yc
 *
 *  glmesh2 functions called:  displayString
 *
 *--------------------------------------------------------------------
 */
static void displayPosition(GLdouble x, GLdouble y)
{
   static GLdouble xsav = 0.0, ysav = 0.0;  /* Previous (x,y) */
   unsigned int fontno = 5;       /* Font for strings */
   GLfloat pos[3];                /* Label position */
   char text[30];                 /* Label text */
   GLdouble ch;                   /* Character height */
   GLboolean xor;

/*
 *  The character height ch is approximately 11 pixels (converted
 *  to object coordinates).  The strings written on the previous
 *  call are erased by writing them in the background color.
 */
   ch = (sx <= sy) ? sx : sy;   /* ch = min(sx, sy) */
   ch = 11.0*ch;
   pos[0] = (GLfloat) (xc-0.5*dx+ch);
   pos[1] = (GLfloat) (yc+0.5*dy-1.5*ch);
   pos[2] = 0.0;
/*
 *  Strings are written in overwrite mode.  The strings written on
 *  the previous call are erased by writing them in the background
 *  color.
 */
   xor = glIsEnabled(GL_COLOR_LOGIC_OP);
   glDisable(GL_COLOR_LOGIC_OP);
   sprintf(text, "x = %.3e", xsav);
   displayString(pos, text, fontno, color_b);
   sprintf(text, "x = %.3e", x);
   displayString(pos, text, fontno, color_t);
   pos[1] = (GLfloat) (yc+0.5*dy-3.0*ch);
   sprintf(text, "y = %.3e", ysav);
   displayString(pos, text, fontno, color_b);
   sprintf(text, "y = %.3e", y);
   displayString(pos, text, fontno, color_t);
/*
 *  Update saved (x,y), and restore the write mode.
 */
   xsav = x;
   ysav = y;
   if (xor) {
      glEnable(GL_COLOR_LOGIC_OP);
      glLogicOp(GL_XOR);
   }
   glutSwapBuffers();
   return;
}
/* End of displayPosition */


/*
 *--------------------------------------------------------------------
 *
 *  displayString:  Write string to the screen with lower left corner
 *                  at raster position pos using font number fontno 
 *                  (defined in the table below) and color scolor.
 *
 *  Note that the raster position is actually a location in the object
 *  coordinate space; i.e., the current modelview and projection 
 *  matrices are applied to pos in order to obtain the screen coordi-
 *  nates of the first character.  If those coordinates are not in the
 *  view volume, nothing is written.  Otherwise the string is written 
 *  but clipped against the window boundaries.
 *
 *  If fontno is invalid, an error message is written to the standard 
 *  output unit and the return value is 1.  Otherwise the return value
 *  is 0.
 *
 *  Global variables required:  dx, dy
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static int displayString(GLfloat pos[3], char *string, 
                         unsigned int fontno, GLfloat scolor[3])
{
   unsigned int i;            /* string index */
   unsigned int nfonts = 9;   /* Number of fonts */
   GLfloat sf;                /* Scale factor for stroked fonts */

   static void *font[9] = {
      GLUT_BITMAP_8_BY_13,         /* 8 x 13 fixed width */
      GLUT_BITMAP_9_BY_15,         /* 9 x 15 fixed width */  
      GLUT_BITMAP_TIMES_ROMAN_10,  /* 10-point proportional spaced */
      GLUT_BITMAP_TIMES_ROMAN_24,  /* 24-point proportional spaced */
      GLUT_BITMAP_HELVETICA_10,    /* 10-point proportional spaced */
      GLUT_BITMAP_HELVETICA_12,    /* 12-point proportional spaced */
      GLUT_BITMAP_HELVETICA_18,    /* 18-point proportional spaced */
      GLUT_STROKE_ROMAN,           /* proportional spaced stroke */
      GLUT_STROKE_MONO_ROMAN       /* mono-spaced Roman simplex */
   };

   if (fontno >= nfonts) {
      printf("displayString:  Invalid font number = %u\n", fontno);
      return 1;
   }
   glColor3fv(scolor);
   if (fontno < 7 ) {
      glRasterPos3fv(pos);
      i = 0;
      while (string[i] != '\0') {
         glutBitmapCharacter(font[fontno], string[i]);
         i = i + 1;
      }
   } else {
      sf = 0.0005*(dx+dy)/4.0;
      glLineWidth(3.0);
      glPushMatrix();
      glTranslatef(pos[0], pos[1], pos[2]);
      glScalef(sf, sf, sf);
      i = 0;
      while (string[i] != '\0') {
         glutStrokeCharacter(font[fontno], string[i]);
         i = i + 1;
      }
      glPopMatrix();
      glLineWidth(line_width);
   }
   return 0;
}
/* End of displayString */


/*
 *--------------------------------------------------------------------
 *
 *  dragCorner:  Mouse callback (callback triggered by a mouse button
 *               press or release with the mouse position in the 
 *               primary window) used to initiate an alteration of the
 *               position of a non-removable vertex.
 *
 *  A left button press may be used to select and drag a non-removable
 *  vertex to any location from which its neighbors are visible.
 *
 *  glmesh2 functions called:  bbox, intsec, xlate
 *
 *--------------------------------------------------------------------
 */
static void dragCorner(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;      /* Computed object coordinates */
   GLboolean intr, visible;
   GLdouble tols, x1, x2, x3, y1, y2, y3;
   GLint i, kt, kv, kv1, kv2, kv3 = 0, kvn, nnb;

   if (button != GLUT_LEFT_BUTTON) return;

   if (state == GLUT_DOWN) {
/*
 *  Compute the object coordinates (wx,wy) of the mouse cursor.
 */
      xlate(x,y,0.0, &wx,&wy,&wz);
/*
 *  Find a non-removable vertex kv1, if any, within squared distance 
 *  tols of the mouse position (wx,wy).
 */
      tols = (sx <= sy) ? sx : sy;   /* tols = min(sx, sy)  */
      tols = 25.0*tols*tols;         /* 5-pixel tolerance   */
      intr = 0;
      for (i = 0; i < nn; i++) {
         kv1 = lnrv[i];
         if ((intr = ((wx-vtxy[kv1][0])*(wx-vtxy[kv1][0]) + 
                      (wy-vtxy[kv1][1])*(wy-vtxy[kv1][1]) <= tols))) 
           break;
      }
      if (!intr) return;
/*
 *  Set up for dragging in Function motion:  store the coordinates of
 *  kv1 in (xv1,yv1), store the indices of the neighboring vertices
 *  of kv1 in the first nnb positions of lste, followed by kv1 in 
 *  position nnb, set XOR write mode, set mode = ALTER_DOMAIN, and
 *  set nnbv = nnb.
 */
      xv1 = vtxy[kv1][0];
      yv1 = vtxy[kv1][1];
      nnb = 0;
      kt = lvtx[kv1];
      do {
         if (ltri[kt][0] == kv1)
            i = 1;
         else if (ltri[kt][1] == kv1)
            i = 2;
         else
            i = 0;
         lste[nnb] = ltri[kt][i];
         nnb++;
         kvn = ltri[kt][(i+1)%3];
         kt = ltri[kt][i+3];
      } while (kt >= 0);
      lste[nnb] = kvn;
      nnb++;
      lste[nnb] = kv1;
      mode = ALTER_DOMAIN;
      nnbv = nnb;
      glEnable(GL_COLOR_LOGIC_OP);
      glLogicOp(GL_XOR);
   } else {
/*
 *  Left button release.
 */
      if (nnbv < 2) return;
      nnb = nnbv;
      kv1 = lste[nnb];
/*
 *  The new point p = (xv1,yv1) is visible iff p Left kv2 -> kv3
 *  for all pairs of adjacent neighbors (kv2,kv3), and there is no
 *  intersection of the path p-kv1 with a boundary edge not containing
 *  kv1 in the same boundary curve.
 */
      visible = 1;
      for (i = 0; i < nnb-1; i++) {
         kv2 = lste[i];
         kv3 = lste[i+1];
         x2 = vtxy[kv2][0];
         y2 = vtxy[kv2][1];
         x3 = vtxy[kv3][0];
         y3 = vtxy[kv3][1];
         if ( (x3-x2)*(yv1-y2) < (y3-y2)*(xv1-x2) ) {
            visible = 0;
            break;   
         }
      }
      if (visible) {
/*
 *  Set kv2 and kv3 to the first and last neighbors of kv1, and 
 *  test for intersections.
 */
         intr = 0;
         kv2 = lste[0];
         x1 = vtxy[kv1][0];
         y1 = vtxy[kv1][1];
         kv = kv2;
         while (kv != kv3) {
            kt = lvtx[kv];
            if (ltri[kt][0] == kv)
               i = 1;
            else if (ltri[kt][1] == kv)
               i = 2;
            else
               i = 0;
            kvn = ltri[kt][i];
            if (intsec(xv1,yv1,x1,y1,vtxy[kv][0],vtxy[kv][1],
                       vtxy[kvn][0],vtxy[kvn][1])) {
               intr = 1;
               break;
            }
            kv = kvn;
         }
         if (intr) visible = 0;
      }
      if (visible) {
/*
 *  Move kv1 to its new position, and update the bounding box.
 */
         vtxy[kv1][0] = xv1;
         vtxy[kv1][1] = yv1;
         bbox();
      }
/*
 *  Set nnbv = 0 indicating that no vertex is currently being dragged.
 */
      nnbv = 0;
      glutPostRedisplay();
   }
   return;
}
/* End of dragCorner */


/*
 *--------------------------------------------------------------------
 *
 *  drawCircle:  Draw a circle of radius cr centered at (cx,cy),
 *               unless cr is too small (less than 5 pixels).
 *
 *  Global variables required:  cr, cx, cy, sx, sy
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void drawCircle(void)
{
   GLdouble c, da, eps, r, s, x, xs, y;
/*
 *  Compute rotation parameters c = cos(da) and s = sin(da) for
 *  angle da = 5/r, where r is the radius in pixels.
 */
   s = (sx >= sy) ? sx : sy;   /* da = max(sx, sy)  */
   r = cr/s;
   if (r < 5.0) return;
   da = 5.0/r;
   c = cos(da);
   s = sin(da);
   eps = cr*(1.0-c)/2.0;
/*
 *  Loop on points uniformly distributed by angle.
 */
   glBegin(GL_LINE_LOOP);
      x = cx + cr;
      y = cy;
      glVertex2d(x, y);
      do {
        xs = x;
        x = cx + c*(x-cx) - s*(y-cy);
        y = cy + s*(xs-cx) + c*(y-cy);
        glVertex2d(x, y);
      } while (x-cx < cr-eps);
   glEnd();
   glutSwapBuffers();
   return;
}
/* End of drawCircle */


/*
 *--------------------------------------------------------------------
 *
 *  drawRect:  Draw the outline of an axis-aligned rectangle in the 
 *             z = 0 plane with lower left corner at (x0,y0), width w, 
 *             and height h.  The origin of the coordinate system is 
 *             at the lower left.
 *
 *  The bottom (lower) and right edges are drawn with color brcolor; 
 *  the top and left edges are drawn with color tlcolor, where brcolor
 *  and tlcolor are arrays of length three containing red, green, and 
 *  blue components normalized to [0,1].  Taking brcolor and tlcolor 
 *  to be black (0,0,0) and white (1,1,1), respectively, gives the
 *  appearance of a button sticking out of the screen.  Reversing 
 *  those colors makes the button appear to extend inward as if it 
 *  were pressed.
 *
 *  The outline is drawn with lines two pixels wide.  If the rectangle
 *  is to be filled, it should be done before the outline is drawn (or
 *  only the interior should be filled) in order to preserve the 
 *  outline colors.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void drawRect(GLuint x0, GLuint y0, GLuint w, GLuint h, 
                     GLfloat brcolor[3], GLfloat tlcolor[3])
{
   glLineWidth(2.0);
   glBegin(GL_LINES);
      glColor3fv(brcolor);
      glVertex2i(x0, y0);            /* Bottom */
      glVertex2i(x0+w-1, y0);

      glVertex2i(x0+w-1, y0);        /* Right */
      glVertex2i(x0+w-1, y0+h-1);

      glColor3fv(tlcolor);
      glVertex2i(x0+w-1, y0+h-1);    /* Top */
      glVertex2i(x0, y0+h-1);

      glVertex2i(x0, y0+h-1);        /* Left */
      glVertex2i(x0, y0);
   glEnd();
   glLineWidth(1.0);
   return;
}
/* End of drawRect */


/*
 *--------------------------------------------------------------------
 *
 *  fillRt:  Fill R with triangles where R contains three adjacent
 *           boundary vertices forming a reentrant angle (corner) at
 *           a removable vertex.
 *
 *  glmesh2 functions called:  bcurve, initR, intsec
 *
 *--------------------------------------------------------------------
 */
static void fillRt(void)
{
   GLboolean intr;
   GLdouble x1, x2, x3, xi, y1, y2, y3, yi;
   GLint i, j, kc, ke, kt, kv1, kv2, kv3, nbv = 0;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
/*
 *  Loop on boundary curves kc.  The sequence of boundary vertex
 *  indices is stored in the first nbv elements of lste.
 */
   for (kc = 0; kc < nc; kc++) {
      bcurve(kc, &nbv, lste);
/*
 *  Inner loop on three consecutive vertices (kv1,kv2,kv3) of 
 *  curve kc.
 */
      for (j = 0; j < nbv; j++) {
         kv1 = lste[j];
         kv2 = lste[(j+1)%nbv];
         kv3 = lste[(j+2)%nbv];
/*
 *  Test for vertices in R.
 */
         if (!listf[kv1]  ||  !listf[kv2]  ||  !listf[kv3]) continue;
/*
 *  Test for kv2 removable.
 */
         intr = 0;
         for (i = 0; i < nn; i++) {
            if (lnrv[i] == kv2) {
               intr = 1;
               break;
            }
         }
         if (intr) continue;
/*
 *  Test for a reentrant corner:  kv2 strictly Left kv1 -> kv3.
 */
         x1 = vtxy[kv1][0];
         y1 = vtxy[kv1][1];
         x2 = vtxy[kv2][0];
         y2 = vtxy[kv2][1];
         x3 = vtxy[kv3][0];
         y3 = vtxy[kv3][1];
         if (!( (x3-x1)*(y2-y1) > (y3-y1)*(x2-x1) )) continue;
/*
 *  Test for no intersection of kv1-kv3 with an edge of the same curve
 *  (that does not contain kv1, kv2, or kv3).
 */
         intr = 0;
         i = (j+3)%nbv;
         while (i != (j-1+nbv)%nbv  &&  nbv > 3) {
            if (intsec(x1,y1,x3,y3,vtxy[lste[i]][0],vtxy[lste[i]][1],
                       vtxy[lste[(i+1)%nbv]][0],
                       vtxy[lste[(i+1)%nbv]][1])) {
               intr = 1;
               break;
            }
            i = (i+1)%nbv;
         }
         if (intr) continue;
         if (nbv > 3) {
/*
 *  A final test for lste[i] not contained in triangle (kv1,kv3,kv2)
 *  is necessary.
 */
            xi = vtxy[lste[i]][0];
            yi = vtxy[lste[i]][1];
            if (( (x3-x1)*(yi-y1) >= (y3-y1)*(xi-x1) )  &&
                ( (x2-x3)*(yi-y3) >= (y2-y3)*(xi-x3) )  &&
                ( (x1-x2)*(yi-y2) >= (y1-y2)*(xi-x2) )) continue;
         }
/*
 *  Add triangle (kv1,kv3,kv2) and edge kv1-kv3 (unless nbv = 3).  
 *  This requires updating ltri, ledg, lvtx, listtr, ntr, lister, 
 *  ner, nt, ne, nb, lcrv, and lste.
 */
         ltri[nt][0] = kv1;
         ltri[nt][1] = kv3;
         ltri[nt][2] = kv2;
         kt = lvtx[kv2];
         if (ltri[kt][0] == kv2)
            i = 2;
         else if (ltri[kt][1] == kv2)
            i = 0;
         else
            i = 1;
         ke = ltri[kt][i+6];
         ledg[ke] = nt;
         ltri[kt][i+3] = nt;
         ltri[nt][3] = kt;
         ltri[nt][6] = ke;

         kt = lvtx[kv1];
         if (ltri[kt][0] == kv1)
            i = 2;
         else if (ltri[kt][1] == kv1)
            i = 0;
         else
            i = 1;
         ke = ltri[kt][i+6];
         ledg[ke] = nt;
         ltri[kt][i+3] = nt;
         ltri[nt][4] = kt;
         ltri[nt][7] = ke;
 
         if (nbv == 3) {
            kt = lvtx[kv3];
            if (ltri[kt][0] == kv3)
               i = 2;
            else if (ltri[kt][1] == kv3)
               i = 0;
            else
               i = 1;
            ke = ltri[kt][i+6];
            ledg[ke] = nt;
            ltri[kt][i+3] = nt;
            ltri[nt][5] = kt;
            ltri[nt][8] = ke;
            listtr[ntr] = nt;
            ntr++;
            nt++;
            nb -= 3;
            for (i = kc; i < nc-1; i++) {
               lcrv[i] = lcrv[i+1];
            }
            nc--;
            break;
         } else {
            ltri[nt][5] = -1;
            ltri[nt][8] = ne;
            ledg[ne] = nt;
            lvtx[kv1] = nt;
            lister[ner] = ne;
            ner++;
            ne++;
            listtr[ntr] = nt;
            ntr++;
            nt++;
            nb--;
            if (lcrv[kc] == kv2) lcrv[kc] = kv1;
            for (i = j+1; i < nbv-1; i++) {   
               lste[i] = lste[i+1];
            }
            nbv--;
            j--;
         }
      }
   }
   return;
}
/* End of fillRt */  


/*
 *--------------------------------------------------------------------
 *
 *  getCenter:  Mouse callback (callback triggered by a mouse button
 *              press or release with the mouse position in the 
 *              primary window) used to move the center of the view
 *              volume.
 *
 *  Set (xc,yc) to the object coordinates corresponding to the mouse 
 *  position when the left button is pressed. 
 * 
 *  glmesh2 functions required:  mouse, reshape, xlate
 *
 *--------------------------------------------------------------------
 */
static void getCenter(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;     /* Computed object coordinates */

   if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
      xlate(x,y,0.0, &wx,&wy,&wz);
      xc = wx;
      yc = wy;
      mode = DEFAULT;
      glutMouseFunc(mouse);            /* Restore standard callback */
      glutSetCursor(GLUT_CURSOR_INHERIT);    /* Restore cursor */
      reshape(glutGet(GLUT_WINDOW_WIDTH),
              glutGet(GLUT_WINDOW_HEIGHT));
      glutPostRedisplay();
   }
   return;
}
/* End of getCenter */


/*
 *--------------------------------------------------------------------
 *
 *  getCircle:  Mouse callback (callback triggered by a mouse button
 *              press or release with the mouse position in the 
 *              primary window) used to select a circle center and
 *              radius.
 *
 *  On left button press, set cr = 0, set (cx,cy) to the object 
 *  coordinates corresponding to the mouse position, set XOR write
 *  mode, and set mode = SELECT_CIRCLE so that Function motion allows
 *  the user to drag a point on the circumference.
 *
 *  On left button release, reset mode = DEFAULT, and restore the 
 *  standard mouse callback, inherited cursor, and overwrite mode.
 *
 *  glmesh2 functions required:  mouse, xlate
 *
 *--------------------------------------------------------------------
 */
static void getCircle(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;     /* Computed object coordinates */

   if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
      xlate(x,y,0.0, &wx,&wy,&wz);
/*
 *  Initialize circle, and set up for dragging in Function motion.
 */
      cr = 0.0;
      cx = wx;
      cy = wy;
      glEnable(GL_COLOR_LOGIC_OP);
      glLogicOp(GL_XOR);
      mode = SELECT_CIRCLE;
   }
   if (button == GLUT_LEFT_BUTTON && state == GLUT_UP) {
/*
 *  Left button release:  reset mode = DEFAULT, and restore the 
 *  default mouse callback, inherited cursor, and overwrite mode.
 */
      mode = DEFAULT;
      glutMouseFunc(mouse);             /* Restore default callback */
      glutSetCursor(GLUT_CURSOR_INHERIT);    /* Restore cursor */
      glDisable(GL_COLOR_LOGIC_OP);
      glutPostRedisplay();
   }
   return;
}
/* End of getCircle */


/*
 *--------------------------------------------------------------------
 *
 *  getCorners:  Replace the set of non-removable vertices lnrv with 
 *               the set of boundary vertices that form corners (do
 *               not lie on the line defined by their predecessors
 *               and successors).
 *
 *  glmesh2 functions called:  bcurve, plDistance
 *
 *--------------------------------------------------------------------
 */
static void getCorners(void)
{
   GLint j, kc, kv1, kv2, kv3, nbv = 0;
   GLdouble tols;
/*
 *  Compute a tolerance (squared distance) for point-line distance.
 */   
   tols = (sx <= sy) ? sx : sy;   /* tols = min(sx, sy)  */
   tols = 25.0*tols*tols;         /* 5-pixel tolerance   */
/*
 *  Initialize nn, and loop on boundary curves kc.  The sequence of
 *  boundary vertex indices is stored in the first nbv elements of 
 *  lste.
 */
   nn = 0;
   for (kc = 0; kc < nc; kc++) {
      bcurve(kc, &nbv, lste);
      for (j = 0; j < nbv; j++) {
         kv1 = lste[(j-1+nbv)%nbv];
	 kv2 = lste[j];
	 kv3 = lste[(j+1)%nbv];
/*
 *  For each subsequence (kv1,kv2,kv3), insert kv2 into lnrv only if
 *  its squared distance from the line kv1-kv3 exceeds tols. 
 */	      
         if ( !(plDistance(vtxy[kv2][0], vtxy[kv2][1], vtxy[kv1][0],
	                   vtxy[kv1][1], vtxy[kv3][0], vtxy[kv3][1],
		           tols) )) {
            lnrv[nn] = kv2;
            nn++;
	 }
      }
   }
   return;
}
/* End of getCorners */  


/*
 *--------------------------------------------------------------------
 *
 *  getLin:  Read a line from stdin into s.
 *
 *  On input:  s is a char array of length at least len.  Input is 
 *             terminated by end of file, a newline character, or when
 *             len-1 characters have been read.
 *
 *  On output:  s contains a (null-terminated) string of length L for 
 *              return value L <= len.  Unlike the C libaray function
 *              fgets, the newline character is not included.  No
 *              characters are left in the input stream.
 *
 *  This function uses the constant EOF defined in <stdio.h>.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static int getLin(char *s, int len)
{
   int c, i;

   i = 0;
   c = getchar();
   while (i < len-1  &&  c != EOF  &&  c != '\n') {
      s[i] = c;
      i++;
      c = getchar();
   }
   s[i] = '\0';
   while (c != EOF && c != '\n') {
      c = getchar();                /* flush input stream */
   }
   return i+1;
}
/* End of getLin */


/*
 *--------------------------------------------------------------------
 *
 *  getPolygon:  Mouse callback (callback triggered by a mouse button
 *               press or release with the mouse position in the 
 *               primary window) used to select a polygonal region R.
 *
 *  On first left button press (np = 0), set (xv1,yv1) and (xv2,yv2)
 *  to the object coordinates corresponding to the mouse position when
 *  the left button is pressed, store (xv1,yv1) in R, and set XOR 
 *  write mode.  On input, mode = CONSTRUCT_R so that the mouse motion 
 *  callback (Function motion) drags a line segment with the mouse 
 *  position.
 *
 *  On subsequent left button presses, set (xv2,yv2) to the mouse
 *  position, and draw the line segment from (xv1,yv1) to (xv2,yv2).
 *
 *  On left button release, set (xv2,yv2) to the mouse position, add
 *  the line segment from (xv1,yv1) to (xv2,yv2) to R, and copy 
 *  (xv2,yv2) to (xv1,yv1) unless the new line segment would result in
 *  a self-intersecting curve.
 *
 *  glmesh2 functions required:  intsec, xlate
 *
 *--------------------------------------------------------------------
 */
static void getPolygon(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;      /* Computed world coordinates */
   GLint i;
   GLboolean intr;

   if (button != GLUT_LEFT_BUTTON) return;
/*
 *  Left button:  compute object coordinates (wx,wy) corresponding
 *                to the mouse position.
 */
   xlate(x,y,0.0, &wx,&wy,&wz);

   if (state == GLUT_DOWN  &&  np == 0) {
/*
 *  First left button press:  store (wx,wy) in (xv1,yv1), (xv2,yv2),
 *  and R.
 */
      xv1 = wx;
      yv1 = wy;
      xv2 = wx;
      yv2 = wy;
      R[np][0] = wx;
      R[np][1] = wy;
      np++;
/*
 *  Set XOR write mode.
 */
      glEnable(GL_COLOR_LOGIC_OP);
      glLogicOp(GL_XOR);
   } 

   if (state == GLUT_DOWN  &&  np > 0) {
/*
 *  Left button press:  set (xv2,yv2) to (wx,wy), and
 *  draw line segment from (xv1,xv2) to (xv2,yv2).
 */
      xv2 = wx;
      yv2 = wy;
      glBegin(GL_LINES);
         glVertex2d(xv1,yv1);
         glVertex2d(xv2,yv2);
      glEnd();
      glutSwapBuffers();
   }

   if (state == GLUT_UP) {
/*
 *  Left button release:  set (xv2,yv2) to (wx,wy).
 */
      xv2 = wx;
      yv2 = wy;
/*
 *  Test for an intersection.
 */
      if (xv2 == xv1  &&  yv2 == yv1) return;
      intr = 0;
      for (i = 0; i <= np-3; i++) {
         if (intsec(xv1,yv1,xv2,yv2,R[i][0],R[i][1],R[i+1][0],
                    R[i+1][1])) {
            intr = 1;
            break;
         }
      }
      if (intr) {
/*
 *  Invalid choice:  erase the line segment, and reset (xv2,yv2) to 
 *  (xv1,yv1).
 */
         glBegin(GL_LINES);
            glVertex2d(xv1,yv1);
            glVertex2d(xv2,yv2);
         glEnd();
         glutSwapBuffers();
         xv2 = xv1;
         yv2 = yv1;
         return;
      }
/*
 *  Store (xv2,yv2) in R, and copy (xv2,yv2) to (xv1,yv1).
 */
      if (np >= npmax) {
         npmax *= 2;
         R = realloc(R, npmax*sizeof *R);
         if (R == NULL) {
            printf("glmesh2:  Unable to allocate sufficient memory"
                   " for R.\n");
            exit(1);
         }
      }
      R[np][0] = xv2;
      R[np][1] = yv2;
      np++;
      xv1 = xv2;
      yv1 = yv2;
   }
   return;
}
/* End of getPolygon */


/*
 *--------------------------------------------------------------------
 *
 *  getTitlePos:  Mouse callback (callback triggered by a mouse button
 *                press or release with the mouse position in the 
 *                primary window) used to select a title position.
 *
 *  This function sets title_pos to the object coordinates of the
 *  mouse position when the left button is pressed. 
 * 
 *  glmesh2 functions required:  mouse, xlate
 *
 *--------------------------------------------------------------------
 */
static void getTitlePos(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;     /* Computed object coordinates */

   if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
      xlate(x,y,0.0, &wx,&wy,&wz);
      title_pos[0] = (GLfloat) wx;
      title_pos[1] = (GLfloat) wy;
      title_pos[2] = (GLfloat) wz;
      mode = DEFAULT;
      glutMouseFunc(mouse);  /* Restore standard callback */
      glutSetCursor(GLUT_CURSOR_INHERIT);  /* Restore cursor */
      glutPostRedisplay();
   }
   return;
}
/* End of getTitlePos */


/*
 *--------------------------------------------------------------------
 *
 *  getVertex:  Mouse callback (callback triggered by a mouse button
 *              press or release with the mouse position in the 
 *              primary window) used to select a vertex in uv-edit
 *              mode.
 *
 *  On left button press, if the mouse position is close enough to a
 *  vertex, write the vertex index, x, y, u, and v values to stdout,
 *  and prompt for new u and v values.
 *
 *  glmesh2 function required:  xlate
 *
 *--------------------------------------------------------------------
 */
static void getVertex(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;      /* Computed world coordinates */
   GLboolean error, intr;
   GLdouble tols;
   GLint kv;
   int  linsiz = 80;
   char line[80];
   char *readerror = "Invalid entry.\n";

   if (button != GLUT_LEFT_BUTTON  ||  state != GLUT_DOWN) return;
/*
 *  Left button press:  compute the object coordinates (wx,wy) of the
 *  mouse position.
 */
   xlate(x,y,0.0, &wx,&wy,&wz);
/*
 *  Find a vertex kv, if any, within squared distance tols of the 
 *  mouse position.
 */
   tols = (sx <= sy) ? sx : sy;   /* tols = min(sx, sy)  */
   tols = 25.0*tols*tols;         /* 5-pixel tolerance   */
   intr = 0;
   for (kv = 0; kv < nv; kv++) {
      if ((intr = ((wx-vtxy[kv][0])*(wx-vtxy[kv][0]) + 
                   (wy-vtxy[kv][1])*(wy-vtxy[kv][1]) <= tols))) 
         break;
   }
   if (!intr) return;
/*
 *  Write kv, x, y, u, and v to stdout, and prompt for new (u,v).
 */
   printf("\n  Vertex %d:  x = %e, y = %e, u = %e, v = %e\n", kv,
          vtxy[kv][0], vtxy[kv][1], uv[kv][0], uv[kv][1]);
   error = 1;
   while (error) {
      printf("  Enter u:  ");
      if (fgets(line, linsiz, stdin) == NULL) return;
      error = (sscanf(line, "%le", &uv[kv][0]) != 1);
      if (error) printf("%s", readerror);
   }
   error = 1;
   while (error) {
      printf("  Enter v:  ");
      if (fgets(line, linsiz, stdin) == NULL) return;
      error = (sscanf(line, "%le", &uv[kv][1]) != 1);
      if (error) printf("%s", readerror);
   }
   printf("  Vertex %d:  x = %e, y = %e, u = %e, v = %e\n\n", kv,
          vtxy[kv][0], vtxy[kv][1], uv[kv][0], uv[kv][1]);
   return;
}
/* End of getVertex */


/*
 *--------------------------------------------------------------------
 *
 *  getViewVolume:  Mouse callback (callback triggered by a mouse 
 *                  button press or release with the mouse position 
 *                  in the primary window) used to select a new view
 *                  volume (pair of corner coordinates).
 *
 *  On left button press, set (xv1,yv1) and (xv2,yv2) to the object
 *  coordinates corresponding to the mouse position when the left 
 *  button is pressed, set XOR write mode, and set mode = ALTER_VIEW_
 *  VOLUME so that the mouse motion callback (Function motion) drags  
 *  a rectangle outline with the mouse position.
 *
 *  On left button release, set (xv2,yv2) to the mouse position, set
 *  (xc,yc), dx, and dy to the center, width, and height of the 
 *  rectangle [xv1,xv2] X [yv1,yv2], restore the standard callback,
 *  cursor, and write mode, and call reshape.
 *
 *  glmesh2 functions required:  mouse, reshape, xlate
 *
 *--------------------------------------------------------------------
 */
static void getViewVolume(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;      /* Computed world coordinates */

   if (button != GLUT_LEFT_BUTTON) return;
/*
 *  Left button:  compute object coordinates (wx,wy) corresponding
 *                to the mouse position.
 */
   xlate(x,y,0.0, &wx,&wy,&wz);

   if (state == GLUT_DOWN) {
/*
 *  Left button press:  set (xv1,yv1) and (xv2,yv2) to (wx,wy).
 */
      xv1 = wx;
      yv1 = wy;
      xv2 = xv1;
      yv2 = yv1;
/*
 *  Set XOR write mode.
 */
      glEnable(GL_COLOR_LOGIC_OP);
      glLogicOp(GL_XOR);
/*
 *  Set mode = ALTER_VIEW_VOLUME and disable the cursor.
 */
      mode = ALTER_VIEW_VOLUME;
      glutSetCursor(GLUT_CURSOR_NONE);
   } else {
/*
 *  Left button release:  reset mode to DEFAULT and set (xv2,yv2) to 
 *  (wx,wy).
 */
      mode = DEFAULT;
      xv2 = wx;
      yv2 = wy;
/*
 *  Restore standard mouse callback, inherited cursor, and overwrite 
 *  mode.
 */
      glutMouseFunc(mouse);
      glutSetCursor(GLUT_CURSOR_INHERIT);
      glDisable(GL_COLOR_LOGIC_OP);
      if (xv2 == xv1  ||  yv2 == yv1) return;
/*
 *  Compute xc, yc, dx, and dy, call reshape, and post redisplay.
 */
      xc = (xv1+xv2)/2.0;
      yc = (yv1+yv2)/2.0;
      dx = xv2-xv1;
      if (dx < 0.0) dx = -dx;
      dy = yv2-yv1;
      if (dy < 0.0) dy = -dy;
      reshape(glutGet(GLUT_WINDOW_WIDTH),glutGet(GLUT_WINDOW_HEIGHT));
      glutPostRedisplay();
   }
   return;
}
/* End of getViewVolume */


/*
 *--------------------------------------------------------------------
 *
 *  init1:  Initialization for window 1.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void init1(void)
{
/*
 *  Set background color.
 */
   glClearColor (color_b[0], color_b[1], color_b[2], 0.0);
   glShadeModel(GL_FLAT);                      /* Default is smooth */
   return;
}
/* End of init1 */


/*
 *--------------------------------------------------------------------
 *
 *  initR:  Initialize user-selected polygon R:  store the lists of
 *          triangles, edges, and vertices contained in R.
 *
 *  Global variables required:
 *
 *      ne, nt, nv = Numbers of edges, triangles, and vertices, 
 *                   respectively, in the mesh.
 *
 *      vtxy = Vertex coordinates.
 *
 *      ltri,ledg = Triangle list and edge list.
 *
 *      newR = Flag with value TRUE iff Function initR needs to be 
 *             called before the triangle mesh can be altered:  set 
 *             to FALSE here.
 *
 *      listf = Array of length nvmax used to store Boolean flags
 *              marking vertices as contained in R.
 *
 *      lister = Array of length 3*nvmax containing the indices of 
 *               the edges that are contained in R (first ner 
 *               positions).
 *
 *      listtr = Array of length 2*nvmax containing the indices of 
 *               the triangles that are contained in R (first ntr 
 *               positions).
 *
 *      listvr = Array of length nvmax containing the indices of the
 *               (triangle mesh) vertices that are contained in R 
 *               (in the first nvr positions).
 *
 *      np = Number of vertices in the boundary of R.
 * 
 *      ner = Number of edges (whose vertices are) contained in R.
 * 
 *      ntr = Number of triangles (whose vertices are) contained in R.
 *
 *      nvr = Number of (triangle mesh) vertices contained in R.
 *
 *      R = Array dimensioned npmax by 2 containing the Cartesian 
 *          coordinates of the sequence of vertices defining the 
 *          boundary of R.
 *
 *  glmesh2 function called:  inside
 *
 *--------------------------------------------------------------------
 */
static void initR(void)
{
   int i, ke, kt, kv, kv1, kv2, kv3;
/*
 *  Test for at least three boundary vertices.
 */
   newR = GL_FALSE;
   if (np < 3) return;
/*
 *  Store nvr and listvr, and mark the vertices that are found in R by
 *  storing flags in listf.
 */
   for (kv = 0; kv < nv; kv++) listf[kv] = 0;
   nvr = 0;
   for (kv = 0; kv < nv; kv++) {
      if (inside(vtxy[kv][0],vtxy[kv][1])) {
         listvr[nvr] = kv;
         nvr++;
         listf[kv] = 1;
      }
   }
/*
 *  Store ntr and listtr, using listf to find the triangles in R.
 */
   ntr = 0;
   for (kt = 0; kt < nt; kt++) {
      kv1 = ltri[kt][0];
      kv2 = ltri[kt][1];
      kv3 = ltri[kt][2];
      if (listf[kv1]  &&  listf[kv2]  &&  listf[kv3]) {
         listtr[ntr] = kt;
         ntr++;
      }
   }
/*
 *  Store ner and lister, using listf to find the edges in R.
 */
   ner = 0;
   for (ke = 0; ke < ne; ke++) {
      kt = ledg[ke];
      if (ltri[kt][6] == ke)
         i = 0;
      else if (ltri[kt][7] == ke)
         i = 1;
      else
         i = 2;
      kv1 = ltri[kt][(i+1)%3];
      kv2 = ltri[kt][(i+2)%3];
      if (listf[kv1]  &&  listf[kv2]) {
         lister[ner] = ke;
         ner++;
      }
   }
   return;
}
/* End of initR */


/*
 *--------------------------------------------------------------------
 *
 *  inputData:  Allocate storage and read a data set containing a tri-
 *              angulation:  number of vertices nv, vertex coordinates
 *              vtxy, number of triangles nt, and CCW-ordered triples
 *              of triangle vertex indices (first three columns of
 *              ltri).
 *
 *  On input:  argc and argv are the parameters input to main.
 *
 *  On output:  Unless the program is terminated with exit code 1
 *              (and an error message), the global variables defining
 *              the triangle mesh are initialized with data from the
 *              input file.  Storage is also allocated for listf, 
 *              lister, listtr, listvr, and lste, and output file
 *              name fname is stored.
 *
 *  Error conditions are tested in the following order:
 *
 *      1:  read error.  This condition is tested following every
 *          read.
 *      2:  nv < 3.
 *      3:  unable to allocate sufficient memory.  This condition
 *          is tested following every call to malloc.
 *      4:  nt < nv-2.
 *      5:  an ltri vertex index is outside its valid range 
 *            (0 to nv-1).
 *      6:  a triangle's vertex indices are not distinct.
 *      7:  two triangles have the same directed edge (pair of 
 *            adjacent vertex indices).
 *      8:  a set of three vertex indices occurs in both clockwise 
 *           and CCW order in ltri.
 *      9:  the same sequence of vertex indices occurs more than once
 *            in ltri (duplicate triangles).
 *      10:  a triangle has three boundary vertices.  This is a 
 *             problem for Dirichlet boundary conditions but not on
 *             an outflow boundary.  Therefore, a warning message is
 *             written, and execution is continued.
 *      11:  an open boundary curve was encountered.
 *      12:  the computed parameters ne, nc, and nb are not consistent
 *             with nv and nt.  This could be caused by a vertex index
 *             appearing more than once in the set of boundary curves.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void inputData(int argc, char *argv[])
{
   int i, i0, i1, j, k, ke, kt, kt0, ktn;
   int  linsiz = 160;
   char line[160];
   FILE *fpr;
   char *filemsg = "\n\n"
     "  glmesh2:  2-D triangle mesh generation using OpenGL.\n\n"
     "  This program requires two command-line arguments:  an input\n"
     "  file name (or path) and an output file name, where the\n"
     "  input file contains the following records:\n\n"
     "    Optional comment lines beginning with # in position 1\n"
     "    nv = number of vertices (at least 3)\n"
     "    x y = Cartesian coordinates of vertex 0\n"
     "    x y = Cartesian coordinates of vertex 1\n"
     "    . . .\n"
     "    x y = Cartesian coordinates of vertex nv-1\n"
     "    Optional comment lines beginning with # in position 1\n"
     "    nt = number of triangles (at least nv-2)\n"
     "    i1 i2 i3 = CCW-ordered vertex indices of triangle 0\n"
     "    i1 i2 i3 = CCW-ordered vertex indices of triangle 1\n"
     "    . . .\n"
     "    i1 i2 i3 = CCW-ordered vertex indices of triangle nt-1\n\n"
     "  Vertex indices are in the range 0 to nv-1.\n";
   char *readerror = "glmesh2:  Error reading file.\n";
   char *memerror = "glmesh2:  Unable to allocate sufficient memory.\n";
/*
 *  Open input file with pointer fpr, and store argv[2] in fname.
 */
   if (argc < 3) {
      printf("%s", filemsg);
      exit(1);
   }
   fpr = fopen(argv[1], "r");
   if (fpr == NULL) {
     printf("glmesh2:  Cannot open file %s.\n", argv[1]);
     exit(1);
   }
   fname = argv[2];
/*
 *  Read a line at a time from fpr into line, beginning with comments.
 */
   line[0] = '#';
   while (line[0] == '#')
      if (fgets(line, linsiz, fpr) == NULL) {
         printf("%s", readerror);
         exit(1);
      }
/*
 *  Read nv. 
 */  
   if (sscanf(line, "%d", &nv) != 1) { 
      printf("%s", readerror);
      exit(1);
   }
   if (nv < 3) {
      printf("glmesh2:  Invalid value:  nv = %d.\n", nv);
      exit(1);
   }
/*
 *  Initialize nvmax to 2*nv, and allocate storage for vtxy, uv, lvtx,
 *  ltri, ledg, lnrv, lcrv, listf, lister, listtr, listvr, and lste.
 */
   nvmax = 2*nv;
   vtxy = malloc (nvmax*sizeof *vtxy);
   uv = malloc (nvmax*sizeof *uv);
   lvtx = malloc (nvmax*sizeof *lvtx);
   ltri = malloc (2*nvmax*sizeof *ltri);
   ledg = malloc (3*nvmax*sizeof *ledg);
   lnrv = malloc (nvmax*sizeof *lnrv);
   lcrv = malloc ((nvmax/3)*sizeof *lcrv);
   listf = malloc (nvmax*sizeof *listf);
   lister = malloc (3*nvmax*sizeof *lister);
   listtr = malloc (2*nvmax*sizeof *listtr);
   listvr = malloc (nvmax*sizeof *listvr);
   lste = malloc (3*nvmax*sizeof *lste);
   if (vtxy == NULL  ||  uv == NULL  ||  lvtx == NULL  ||  
       ltri == NULL  ||  ledg == NULL  ||  lnrv == NULL  ||  
       lcrv == NULL  ||  listf == NULL  ||  lister == NULL  ||  
       listtr == NULL  ||  listvr == NULL  ||  lste == NULL) {
      printf("%s", memerror);
      exit(1);
   }
/*
 *  Read vertices vtxy.
 */
   for (i = 0; i < nv; i++) {
      if (fgets(line, linsiz, fpr) == NULL) {
         printf("%s", readerror);
         exit(1);
      }
      if (sscanf(line, "%lf%lf", &vtxy[i][0], &vtxy[i][1]) != 2) {
         printf("%s", readerror);
         exit(1);
      }
   }
/*
 *  Store zeros in uv.
 */
   for (i = 0; i < nv; i++) {
      uv[i][0] = 0.0;
      uv[i][1] = 0.0;
   }
/*
 *  Read comment lines if any.
 */
   line[0] = '#';
   while (line[0] == '#')
      if (fgets(line, linsiz, fpr) == NULL) {
         printf("%s", readerror);
         exit(1);
      }
/*
 *  Read number of triangles nt.
 */
   if (sscanf(line, "%d", &nt) != 1) {
      printf("%s", readerror);
      exit(1);
   }
   if (nt < nv-2) {
      printf("glmesh2:  Invalid value:  nt = %d.\n", nt);
      exit(1);
   }
/*
 *  Read triangle list ltri.
 */
   for (kt = 0; kt < nt; kt++) {
      if (fgets(line, linsiz, fpr) == NULL) {
         printf("%s", readerror);
         exit(1);
      }
      if (sscanf(line, "%u%u%u", &ltri[kt][0], &ltri[kt][1],
                                  &ltri[kt][2]) != 3) {
         printf("%s", readerror);
         exit(1);
      }
      if (ltri[kt][0] < 0 || ltri[kt][0] >= nv ||
          ltri[kt][1] < 0 || ltri[kt][1] >= nv ||
          ltri[kt][2] < 0 || ltri[kt][2] >= nv) {
         printf("glmesh2:  The i-th triple includes an invalid index "
                "for i = %d.\n", kt);
         exit(1);
      }
      if (ltri[kt][0] == ltri[kt][1] || ltri[kt][1] == ltri[kt][2] ||
          ltri[kt][2] == ltri[kt][0]) {
         printf("glmesh2:  The vertex indices are not distinct in the "
                "i-th triple for i = %d.\n", kt);
         exit(1);
      }
   }
/*
 *  Initialize ltri columns 3 to 8 with flag values -1 to indicate
 *  values that have not yet been computed, and initialize counts.
 */
   for (kt = 0; kt < nt; kt++) {
      for (j = 3; j < 9; j++) {
         ltri[kt][j] = -1; 
      }
   }
   nc = 0;
   nb = 0;
   ne = 0;
/*
 *  Construct the remaining portion of ltri by searching for
 *  neighboring triangles and generating edge indices.  The outer 
 *  loop is on triangles kt;  inner loop on potential neighbors ktn.
 */
   for (kt = 1; kt < nt; kt++) {
      for (ktn = 0; ktn < kt; ktn++) {
/*
 *  Loop on pairs of edges ltri[kt][i]-ltri[kt][i+1],
 *  ltri[ktn][j]-ltri[ktn][j+1].
 */
         for (i = 0; i < 3; i++) {
            for (j = 0; j < 3; j++) {
/*
 *  Test for error 7.
 */
               if (ltri[kt][i] == ltri[ktn][j]  &&
                   ltri[kt][(i+1)%3] == ltri[ktn][(j+1)%3]) {
                  printf("glmesh2:  Two or more triangles have the "
                         "same pair of adjacent vertex indices.\n");
                  exit(1);
               }
               if (ltri[kt][i] != ltri[ktn][(j+1)%3]  ||
                   ltri[kt][(i+1)%3] != ltri[ktn][j]) continue;
/*
 *  Shared edge:  ltri[kt][i]-ltri[kt][i+1] =
 *                ltri[ktn][j+1]-ltri[ktn][j]
 *
 *  Test for errors 7, 8, and 9, and store the indices.
 */
               if (ltri[kt][(i+2)%3] == ltri[ktn][(j+2)%3]) {
                  printf("glmesh2:  The same triple of vertex indices "
                         "occurs\n with both cyclical orderings.\n");
                  exit(1);
               }
               if (ltri[ktn][(j+2)%3+3] >= 0) {
                  k = ltri[ktn][(j+2)%3+3];
                  if (ltri[kt][(i+2)%3] == ltri[k][(i+2)%3]) {
                     printf("glmesh2:  Duplicate vertex index "
                            "triples.\n");
                     exit(1);
                  }
                  printf("glmesh2:  Two or more triangles have the "
                         "same pair of adjacent vertex indices.\n");
                  exit(1);
               }
               ltri[kt][(i+2)%3+3] = ktn;
               ltri[ktn][(j+2)%3+3] = kt;
               ne++;
               ltri[kt][(i+2)%3+6] = ne-1;
               ltri[ktn][(j+2)%3+6] = ne-1;
            }
         }
      }
   }
/*
 *  Count the boundary curves nc and boundary edges (vertices) nb, 
 *  assign them edge indices, and construct lcrv.
 */
   for (kt0 = 0; kt0 < nt; kt0++) {
      for (j = 0; j < 3; j++) {
         if (ltri[kt0][j+6] >= 0) continue;
/*
 *  New boundary curve beginning with vertices i0, i1.
 *  Initialize for traversal.  A counter k is used to
 *  avoid an infinite loop.
 */
         kt = kt0;
         i = j;
         nc++;
         nb++;
         ne++;
         ltri[kt][i+6] = ne-1;
         i = (i+1)%3;
         i1 = ltri[kt][i];
         i = (i+1)%3;
         i0 = ltri[kt][i];
         lcrv[nc-1] = i0;
         k = 1;
/*
 *  Boundary curve traversal:  find the next vertex.
 */
         while (i1 != i0) {
            ktn = ltri[kt][i+3];
            while (ktn >= 0) {
               kt = ktn;
               if (ltri[kt][0] == i1) 
                  i = 1;
               else if (ltri[kt][1] == i1) 
                  i = 2;
               else 
                  i = 0;
               ktn = ltri[kt][i+3];
            }
            if (kt == kt0) 
               printf("glmesh2:  WARNING:  Triangle %d has three "
                      "boundary vertices.\n", kt);
            kt0 = kt;
            nb++;
            ne++;
            ltri[kt][i+6] = ne-1;
            i = (i+1)%3;
            i1 = ltri[kt][i];
            i = (i+1)%3;
            k++;
            if (k > nv) {
               printf("glmesh2:  Open boundary curve.\n");
               exit(1);
            }
         }
      }
   }
/*
 *  Test for error 12.
 */
   if (nb < 3*nc  ||
       ne != 3*nv-nb+3*nc-6  ||
       nt != 2*nv-nb+2*nc-4) {
      printf("glmesh2:  Inconsistent counts:  nv = %d, nb = %d, "
             "nc = %d, nt = %d, ne = %d\n", nv, nb, nc, nt, ne);
      exit(1);
   }
/*
 *  Construct lvtx.  Initialize the list with flag values -1
 *  indicating entries not yet set.
 */
   for (i0 = 0; i0 < nv; i0++)
      lvtx[i0] = -1;
/*
 *  Store triangle indices of boundary vertices i0 in a loop 
 *  on triangle/vertex pairs.
 */
   for (kt = 0; kt < nt; kt++) {
      for (i = 0; i < 3; i++) {
         if (ltri[kt][i+3] < 0) {
            i0 = ltri[kt][(i+1)%3];
            lvtx[i0] = kt;
         }
      }
   }
/*
 *  Fill in triangle indices for vertices with negative entries,
 *  and store ledg.
 */
   for (kt = 0; kt < nt; kt++) {
      for (i = 0; i < 3; i++) {
         i0 = ltri[kt][i];
         if (lvtx[i0] < 0) lvtx[i0] = kt;
         if (ltri[kt][i+3] < kt) {
            ke = ltri[kt][i+6];
            ledg[ke] = kt;
         }
      }
   }
   return;
} 
/* End of inputData */


/*
 *--------------------------------------------------------------------
 *
 *  insertb:  Partition a triangle into three equal-area subtriangles
 *            by inserting the barycenter as a new vertex.  The uv
 *            values at the barycenter are computed by linear
 *            interpolation.
 *
 *  On input:  kt0 is the index of the triangle to be partitioned.
 *
 *  On output:  Unless nv >= nvmax or kt0 is an invalid index or kt0
 *              does not index a triangle in R, the triangle mesh and
 *              polygon lists for R are updated with the revision.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void insertb(int kt0)
{
   int i, ke1, ke2, ke3, kt1, kt2, kt3, kv1, kv2, kv3;

/*
 *  Test for invalid input.
 */
   if (nv >= nvmax  ||  kt0 < 0  ||  kt0 >= nt) return;
/*
 *  Store indices of vertices, opposite triangles, and edges of kt0
 *  in local variables.
 */
   kv1 = ltri[kt0][0];
   kv2 = ltri[kt0][1];
   kv3 = ltri[kt0][2];
   if (!listf[kv1]  ||  !listf[kv2]  ||  !listf[kv3]) return;
   kt1 = ltri[kt0][3];
   kt2 = ltri[kt0][4];
   kt3 = ltri[kt0][5];
   ke1 = ltri[kt0][6];
   ke2 = ltri[kt0][7];
   ke3 = ltri[kt0][8];
/*
 *  Compute the barycenter with index nv, and update vtxy and uv.
 */
   vtxy[nv][0] = (vtxy[kv1][0] + vtxy[kv2][0] + vtxy[kv3][0])/3.0;
   vtxy[nv][1] = (vtxy[kv1][1] + vtxy[kv2][1] + vtxy[kv3][1])/3.0;
   uv[nv][0] = (uv[kv1][0] + uv[kv2][0] + uv[kv3][0])/3.0;
   uv[nv][1] = (uv[kv1][1] + uv[kv2][1] + uv[kv3][1])/3.0;
/*
 *  Overwrite row kt0, append two new rows to ltri, and adjust
 *  additional ltri entries.
 */
   ltri[kt0][0] = nv;
   ltri[kt0][4] = nt;
   ltri[kt0][5] = nt+1;
   ltri[kt0][7] = ne+2;
   ltri[kt0][8] = ne+1;

   ltri[nt][0] = nv;
   ltri[nt][1] = kv3;
   ltri[nt][2] = kv1;
   ltri[nt][3] = kt2;
   ltri[nt][4] = nt+1;
   ltri[nt][5] = kt0;
   ltri[nt][6] = ke2;
   ltri[nt][7] = ne;
   ltri[nt][8] = ne+2;

   ltri[nt+1][0] = nv;
   ltri[nt+1][1] = kv1;
   ltri[nt+1][2] = kv2;
   ltri[nt+1][3] = kt3;
   ltri[nt+1][4] = kt0;
   ltri[nt+1][5] = nt;
   ltri[nt+1][6] = ke3;
   ltri[nt+1][7] = ne+1;
   ltri[nt+1][8] = ne;
   if (kt2 >= 0) {
      if (ltri[kt2][3] == kt0)
         i = 3;
      else if (ltri[kt2][4] == kt0)
         i = 4;
      else
         i = 5;
      ltri[kt2][i] = nt;
   }
   if (kt3 >= 0) {
      if (ltri[kt3][3] == kt0)
         i = 3;
      else if (ltri[kt3][4] == kt0)
         i = 4;
      else
         i = 5;
      ltri[kt3][i] = nt+1;
   }
/*
 *  Adjust ledg and lvtx entries.
 */
   ledg[ke2] = nt;
   ledg[ke3] = nt+1;
   ledg[ne] = nt+1;
   ledg[ne+1] = nt+1;
   ledg[ne+2] = nt;
   if (lvtx[kv1] == kt0) lvtx[kv1] = nt+1;
   if (lvtx[kv3] == kt0) lvtx[kv3] = nt;
   lvtx[nv] = kt0;
/*
 *  Append new elements to lists associated with R.
 */
   listvr[nvr] = nv;
   nvr++;
   lister[ner] = ne;         
   lister[ner+1] = ne+1;
   lister[ner+2] = ne+2;
   ner += 3;
   listtr[ntr] = nt;
   listtr[ntr+1] = nt+1;
   ntr += 2;
   listf[nv] = 1;
/*
 *  Update vertex, edge, and triangle counts.
 */
   nv++;
   ne += 3;
   nt += 2;
   return;
}
/* End of insertb */


/*
 *--------------------------------------------------------------------
 *
 *  inside:  Locate a point P relative to the polygon R, returning a
 *           nonzero value iff P is contained in R.
 *
 *  R is defined by a cyclically ordered sequence of vertices defining
 *  a simple closed curve.  Adjacent vertices need not be distinct, 
 *  and the curve orientation is irrelevant. 
 *
 *  The algorithm is referred to as ray casting, crossing number, or
 *  even-odd rule.  The point P is inside R iff there is an odd number
 *  of intersections between the polygon boundary curve and a hori-
 *  zontal ray with endpoint P.  Intersections are counted by testing
 *  each curve segment (polygon side) for intersection.  In order to 
 *  avoid counting a vertex as two intersections when it lies on the 
 *  ray and has neighbors on opposite sides of the ray, we arbitrarily
 *  define vertices with the same y component as P as being above the
 *  ray.
 *
 *  If P lies on the boundary of R, the decision is unspecified, and
 *  if P is close to the polygon boundary, the problem is ill-
 *  conditioned, and the decision may be incorrect. 
 *
 *  On input:  xp,yp = Cartesian coordinates of P.
 *
 *  Global variables required:
 *
 *      np = Number of vertices in the boundary of R.
 *
 *      R = Array dimensioned np by 2 containing the Cartesian 
 *          coordinates of the sequence of vertices defining the 
 *          boundary of R.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static unsigned char inside(double xp, double yp)
{
   unsigned char odd = 0;
   int i;
   double x, x1, x2, y1, y2;

   x2 = R[np-1][0];
   y2 = R[np-1][1];
   for (i = 0; i < np; i++) {
      x1 = x2;
      y1 = y2;
      x2 = R[i][0];
      y2 = R[i][1];
/*
 *  Test for an intersection with the line y = yp.
 */
      if ((y1 < yp  &&  y2 < yp)  ||
          (y1 >= yp  &&  y2 >= yp)) continue;
/*
 *  Compute the x component of the point of intersection, and test
 *  for an intersection with the ray x > xp.
 */
      x = x1 + (yp-y1)*(x2-x1)/(y2-y1);
      if (x > xp) odd = !odd;
   }
   return odd;
}
/* End of inside */


/*
 *--------------------------------------------------------------------
 *
 *  intsec:  Test for an intersection between a pair of line 
 *           segments in the plane.
 *
 *  Note that an incorrect decision may result from floating point 
 *  error if the four endpoints are nearly collinear.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static unsigned char intsec(double x1, double y1, double x2, 
                            double y2, double x3, double y3, 
                            double x4, double y4)
{
   double a, b, d, dx12, dx31, dx34, dy12, dy31, dy34;
/*
 *  Test for overlap between the smallest rectangles that contain
 *  the line segments and have sides parallel to the axes.
 */
   if ((x1 < x3  &&  x1 < x4  &&  x2 < x3  &&  x2 < x4)  ||
       (x1 > x3  &&  x1 > x4  &&  x2 > x3  &&  x2 > x4)  ||
       (y1 < y3  &&  y1 < y4  &&  y2 < y3  &&  y2 < y4)  ||
       (y1 > y3  &&  y1 > y4  &&  y2 > y3  &&  y2 > y4)) {
      return 0;
   }
/*
 *  Compute a = (p4-p3) X (p1-p3), b = (p2-p1) X (p1-p3), and
 *  d = (p2-p1) X (p4-p3) (z components).
 */
   dx12 = x2 - x1;
   dy12 = y2 - y1;
   dx34 = x4 - x3;
   dy34 = y4 - y3;
   dx31 = x1 - x3;
   dy31 = y1 - y3;
   a = dx34*dy31 - dx31*dy34;
   b = dx12*dy31 - dx31*dy12;
   d = dx12*dy34 - dx34*dy12;
/*
 *  If d == 0, either the line segments are parallel, or one (or both)
 *  of them is a single point.
 */
   if (d == 0.0) return (a == 0.0  &&  b == 0.0);
/*
 *  d != 0 and the point of intersection of the lines defined by the 
 *  line segments is p = p1 + (a/d)*(p2-p1) = p3 + (b/d)*(p4-p3).
 */
   return (a/d >= 0.0  &&  a/d <= 1.0  &&
           b/d >= 0.0  &&  b/d <= 1.0);   
}
/* End of intsec */


/*
 *--------------------------------------------------------------------
 *
 *  key:  Associate ASCII keys with menu entries.
 *        The mouse coordinates (x,y) are not used.
 *
 *  glmesh2 function called:  menu
 *
 *--------------------------------------------------------------------
 */
static void key(unsigned char key, int x, int y)
{
   menu((int) key);
   return;
}
/* End of key */


/*
 *--------------------------------------------------------------------
 *
 *  makeMenu:  Create a menu (for window 1) attached to the right 
 *             mouse button.
 *
 *  glmesh2 function required:  menu
 *
 *--------------------------------------------------------------------
 */
static void makeMenu(void)
{
   int submenu1, submenu2;

   glutCreateMenu(menu);

   glutAddMenuEntry("", 0);
   glutAddMenuEntry(" a:  Toggle preservation of Aspect ratio", 'a');
   glutAddMenuEntry(" b:  Toggle display of Bounding box", 'b');
   glutAddMenuEntry(" B:  Set R to contain the Bounding box", 'B');
   glutAddMenuEntry(" c:  Toggle triangle color-coding by shape quality", 'c');
   glutAddMenuEntry(" e:  Toggle annotation of Edges with indices", 'e');
   glutAddMenuEntry(" F:  Output the triangle mesh to a File", 'F');
   glutAddMenuEntry(" i:  Triangulation Integrity check (debugging)", 'i');
   glutAddMenuEntry(" j:  Toggle antialiasing (Jaggies)", 'j');
   glutAddMenuEntry(" K:  Increase line width", 'K');
   glutAddMenuEntry(" k:  Decrease line width", 'k');
   glutAddMenuEntry(" o:  Draw a circle (drag with left mouse button)", 'o');
   glutAddMenuEntry(" p:  Pan (select new center with mouse button)", 'p');
   glutAddMenuEntry(" P:  Print triangle mesh data structure", 'P');
   glutAddMenuEntry(" q:  Query current state (print parameters)", 'q');
   glutAddMenuEntry(" r:  Restore default view volume", 'r');
   glutAddMenuEntry(" R:  Select a new polygon R (mouse button presses)", 'R');
   glutAddMenuEntry(" t:  Toggle annotation of Triangles with indices", 't');
   glutAddMenuEntry(" u:  Enter uv-edit mode (select vertex with mouse)", 'u');
   glutAddMenuEntry(" v:  Toggle annotation of Vertices with indices", 'v');
   glutAddMenuEntry(" V:  Toggle annotation of non-removable Vertices", 'V');
   glutAddMenuEntry(" w:  Select a new Window (left mouse button)", 'w');
   glutAddMenuEntry(" z:  Zoom in", 'z');
   glutAddMenuEntry(" Z:  Zoom out", 'Z');
   glutAddMenuEntry(" [Arrows]:  Move center", 0);
   glutAddMenuEntry(" [Home]:  Restore default view volume", 0);
   glutAddMenuEntry(" [Page Up, Page Down]:  Zoom in, zoom out", 0);
   glutAddMenuEntry(" [Left mouse button]:  Alter polygon R", 0);
   glutAddMenuEntry(" [Shifted left mouse button]:  Delete vertex of R", 0);
   glutAddMenuEntry(" [Enter]:  Terminate mode (D, E, or R)", 13);
   glutAddMenuEntry(" [esc]:  Cancel mode (D/E/R), or terminate program", 27);

   submenu1 = glutCreateMenu(menu);
   glutAddMenuEntry(" C:  Replace non-removable vertex set by corners", 'C');
   glutAddMenuEntry(" D:  Alter Domain (drag non-removable vertices)", 'D');
   glutAddMenuEntry(" E:  Select Edges for swapping (left mouse button)", 'E');
   glutAddMenuEntry(" f:  Fill R with triangles", 'f');
   glutAddMenuEntry(" L:  Toggle boundary Layer (smoothing method)", 'L');
   glutAddMenuEntry(" n:  Convert vertices in R to removable", 'n');
   glutAddMenuEntry(" N:  Convert vertices in R to Non-removable", 'N');
   glutAddMenuEntry(" O:  Move non-removable vertices in R to circle", 'O');
   glutAddMenuEntry(" s:  Optimize R with Delaunay edge Swaps", 's');
   glutAddMenuEntry(" S:  Optimize R by Smoothing (moving vertices)", 'S');
   glutAddMenuEntry(" +:  Refine R by adding triangle barycenters", '+');
   glutAddMenuEntry(" >:  Refine R by adding edge midpoints", '>');
   glutAddMenuEntry(" -:  Unrefine by reversing + or > if possible", '-');
   glutAddMenuEntry(" <:  Coarsen R by removing interior vertices", '<');
   glutAddMenuEntry(" [Backspace or Delete]:  Remove triangles in R", 8);

   submenu2 = glutCreateMenu(menu);
   glutAddMenuEntry(" T:  Replace Text (keyboard string entry)", 'T');
   glutAddMenuEntry(" l:  Change title Location (left button)", 'l');
   glutAddMenuEntry(" !:  Toggle font selection:  indices/title", '!');
   glutAddMenuEntry(" 0:  Set 8 by 13 fixed width font", '0');
   glutAddMenuEntry(" 1:  Set 9 by 15 fixed width font", '1');
   glutAddMenuEntry(" 2:  Set 10-point prop. spaced Times Roman", '2');
   glutAddMenuEntry(" 3:  Set 24-point prop. spaced Times Roman", '3');
   glutAddMenuEntry(" 4:  Set 10-point prop. spaced Helvetica", '4');
   glutAddMenuEntry(" 5:  Set 12-point prop. spaced Helvetica", '5');
   glutAddMenuEntry(" 6:  Set 18-point prop. spaced Helvetica", '6');
   glutAddMenuEntry(" 7:  Set prop. spaced Roman stroke font", '7');
   glutAddMenuEntry(" 8:  Set mono-spaced Roman stroke font", '8');

   glutSetMenu(1);
   glutAddSubMenu(" Triangulation Updates -->", submenu1);
   glutAddSubMenu(" Title and Font -->", submenu2);

   glutAttachMenu(GLUT_RIGHT_BUTTON);
   return;
}
/* End of makeMenu */


/*
 *--------------------------------------------------------------------
 *
 *  menu:  Respond to menu selections.
 *
 *  glmesh2 functions required:  all of them
 *
 *--------------------------------------------------------------------
 */
static void menu(int item)
{
   static unsigned char font_select = 0;

   switch (item) {
   case 'a':   /* Toggle preservation of aspect ratio */
      aspect = !aspect;
      reshape(glutGet(GLUT_WINDOW_WIDTH),
              glutGet(GLUT_WINDOW_HEIGHT));
      break;

   case 'B':   /* Set R to a rectangle containing the bounding box */
      np = 4;
      R[0][0] = xmin - 0.05*(xmax-xmin);
      R[0][1] = ymin - 0.05*(ymax-ymin);
      R[1][0] = xmax + 0.05*(xmax-xmin);
      R[1][1] = ymin - 0.05*(ymax-ymin);
      R[2][0] = xmax + 0.05*(xmax-xmin);
      R[2][1] = ymax + 0.05*(ymax-ymin);
      R[3][0] = xmin - 0.05*(xmax-xmin);
      R[3][1] = ymax + 0.05*(ymax-ymin);
      initR();
      break;

   case 'b':   /* Toggle display of bounding box */
      plot_box = !plot_box;
      break;

   case 'c':   /* Toggle color-coding by a triangle quality measure */
      color_code = !color_code;
      break;
      
   case 'e':   /* Toggle display of edge indices */
      plot_e = !plot_e;
      break;

   case 'F':   /* Write the triangle mesh to the output file */
      trwrite(fname, nc, nb, nv, nt, ne, vtxy, uv, ltri);
      break;

   case 'i':   /* Triangulation integrity check */
      trmtst(nc, nb, nv, nn, nt, ne, ltri, lvtx, ledg, lnrv, lcrv);
      return;

   case 'j':   /* Toggle antialiasing */
      antialias = !antialias;
      if (antialias) {
         glHint(GL_LINE_SMOOTH_HINT, GL_NICEST);
         glEnable(GL_LINE_SMOOTH);
         glEnable(GL_BLEND);
         glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
      } else {
         glDisable(GL_LINE_SMOOTH);
         glDisable(GL_BLEND);
      }
      break;

   case 'K':   /* Increase line width */
      line_width *= 1.25;
      break;

   case 'k':   /* Decrease line width */
      line_width *= 0.80;
      break;

   case 'o':   /* Draw a circle by dragging a point on the 
                  circumference with the left mouse button */
      glutMouseFunc(getCircle);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      break;

   case 'p':   /* Pan:  select center (xc,yc) */
      mode = PAN;
      glutMouseFunc(getCenter);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      break;

   case 'P':   /* Print triangle mesh data structure */
      trprint(nc, nb, nv, nn, nt, ne, vtxy, uv, ltri, lvtx, ledg, 
              lnrv, lcrv);
      return;

   case 'q':   /* Query current state */
      printf("\n\n");
      printf("  Mesh parameters:  nc = %u, nb = %u, nv = %u, nn = %u"
             ", nt = %u, ne = %u\n", nc, nb, nv, nn, nt, ne);
      printf("                    nvmax = %u\n\n", nvmax);
      printf("  View volume parameters:\n"
             "     Center x = %g\n", xc);
      printf("     Center y = %g\n", yc);
      printf("     Width = %g\n", dx);
      printf("     Height = %g\n\n", dy);
      printf("  Boundary layer flag = %s\n", flagName[bdry_layer]);
      printf("  Antialiasing flag = %s\n", flagName[antialias]);
      printf("  Preserve aspect ratio flag = %s\n", flagName[aspect]);
      printf("  Line width = %.3f\n", line_width);
      printf("  Font number for indices = %u\n", ind_font);
      printf("  Font number for title = %u\n", title_font);
      printf("  Number of vertices in boundary of R = %d\n", np);
      printf("  Current mode = %s\n", modeName[mode]);
      return;

   case 'r':   /* Restore default view volume */
      dx = 1.2*(xmax-xmin);
      dy = 1.2*(ymax-ymin);
      xc = (xmin+xmax)/2.0;
      yc = (ymin+ymax)/2.0;
      reshape(glutGet(GLUT_WINDOW_WIDTH),
              glutGet(GLUT_WINDOW_HEIGHT));
      break;
    
   case 'R':   /* Select a new polygon R */
      mode = CONSTRUCT_R;
      newR = GL_TRUE;
      np = 0;
      ner = 0;
      ntr = 0;
      nvr = 0;
      glutMouseFunc(getPolygon);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      break;

   case 't':   /* Toggle display of triangle indices */
      plot_t = !plot_t;
      break;
    
   case 'u':   /* Enter uv-edit mode */
      if (mode == CONSTRUCT_R) {   /* Terminate polygon selection  */
         glDisable(GL_COLOR_LOGIC_OP);
         np = 0;
         ner = 0;
         ntr = 0;
         nvr = 0;
      }
      if (mode == ALTER_DOMAIN) {  /* Terminate corner drag mode */
         glDisable(GL_COLOR_LOGIC_OP);
      }
      mode = EDIT_UV;
      glutMouseFunc(getVertex);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      return;

   case 'v':   /* Toggle display of vertex indices */
      plot_v = !plot_v;
      break;
    
   case 'V':   /* Toggle display of non-removable vertex indices */
      plot_nrv = !plot_nrv;
      break;

   case 'w':   /* Select view volume center, width, and height */
      glutMouseFunc(getViewVolume);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      return;

   case 'z':   /* Zoom in */
      dx *= 0.8;
      dy *= 0.8;
      reshape(glutGet(GLUT_WINDOW_WIDTH),
              glutGet(GLUT_WINDOW_HEIGHT));
      break;

   case 'Z':   /* Zoom out */
      dx *= 1.25;
      dy *= 1.25;
      reshape(glutGet(GLUT_WINDOW_WIDTH),
              glutGet(GLUT_WINDOW_HEIGHT));
      break;

   case 10:    /* Terminate current mode */
   case 13: 
      if (mode == CONSTRUCT_R) closeR();
      mode = DEFAULT;          /* Restore default mode */
      glutMouseFunc(mouse);    /* Restore default mouse callback */
      glutSetCursor(GLUT_CURSOR_INHERIT);  /* Standard cursor */
      glDisable(GL_COLOR_LOGIC_OP);        /* Overwrite mode */
      break;
      
   case 27:    /* Terminate program or current mode */
      if (mode == CONSTRUCT_R) {
         mode = DEFAULT;        /* Terminate polygon selection mode */
         glutMouseFunc(mouse);
         glutSetCursor(GLUT_CURSOR_INHERIT);
         glDisable(GL_COLOR_LOGIC_OP);
         np = 0;
         ner = 0;
         ntr = 0;
         nvr = 0;
         break;
      }
      if (mode == ALTER_DOMAIN  ||  mode == SWAP_EDGES  ||  
          mode == EDIT_UV) {
         mode = DEFAULT;        /* Terminate drag mode */
         glutMouseFunc(mouse);
         glutSetCursor(GLUT_CURSOR_INHERIT);
         glDisable(GL_COLOR_LOGIC_OP);
         break;
      }
      sprintf(dlg_string,"You really want to quit?");
      dlg_exit = GL_TRUE;
      glutHideWindow();          /* Hide window 1 */
      glutSetWindow(window[1]);  /* Make window 2 current */
      glutShowWindow();          /* Make window 2 visible */
      dlgDisplay();
      break;

/*
 *  Submenu 1:
 */
   case 'C':   /* Replace set of non-removable vertices lnrv by 
                  the set of corners */
      getCorners();
      break;

   case 'D':   /* Drag a non-removable vertex */
      if (mode == CONSTRUCT_R) {   /* Terminate polygon selection */
         np = 0;
         ner = 0;
         ntr = 0;
         nvr = 0;
      }
      if (newR) initR();
      mode = ALTER_DOMAIN;
      glutMouseFunc(dragCorner);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      break;

   case 'E':   /* Enter edge-swapping mode (left button selection) */
      if (mode == CONSTRUCT_R) {   /* Terminate polygon selection  */
         glDisable(GL_COLOR_LOGIC_OP);
         np = 0;
         ner = 0;
         ntr = 0;
         nvr = 0;
      }
      if (mode == ALTER_DOMAIN) {  /* Terminate corner drag mode */
         glDisable(GL_COLOR_LOGIC_OP);
      }
      mode = SWAP_EDGES;
      glutMouseFunc(swapEdges);
      glutSetCursor(GLUT_CURSOR_CROSSHAIR);  /* Change cursor */
      return;

   case 'f':   /* Fill R with triangles */
      fillRt();
      break;

   case 'L':   /* Toggle boundary layer smoothing */
      bdry_layer = !bdry_layer;
      return;

   case 'n':   /* Convert vertices in R to removable */
      convertRv(1);
      break;

   case 'N':   /* Convert vertices in R to non-removable */
      convertRv(0);
      break;

   case 'O':   /* Project non-removable vertices in R onto circle */
      moveRn();
      break;

   case 's':   /* Optimize R with edge swaps */
      optimizeRt();
      break;

   case 'S':   /* Optimize R by smoothing (moving vertices) */
      optimizeRv();
      break;

   case '+':   /* Refine R by adding triangle barycenters */
      refineRt();
      break;

   case '>':   /* Refine R by adding edge midpoints */
      refineRe();
      break;

   case '-':   /* Unrefine the mesh by reversing the previous 
                  refinement by refineRe or refineRt */
      unrefine();
      break;

   case '<':   /* Coarsen R by removing interior vertices */
      coarsenRv();
      break;

   case 8:     /* Remove triangles in R */
   case 127:   /* Remove triangles in R */
      removeRt();
      break;

/*
 *  Submenu 2:
 */
   case 'T':   /* Input a new character string for title */
      printf("\n  Enter a title string:  ");
      getLin(title, LSTRING);
      break;

   case 'l':   /* Get a new title location title_pos */
      mode = PLACE_TITLE;
      glutMouseFunc(getTitlePos);
      glutSetCursor(GLUT_CURSOR_RIGHT_ARROW);  /* Change cursor */
      break;
 
   case '!':   /* Toggle font selection ('0' to '8') between indices
                  (font_select = 0) and title (font_select = 1) */
      font_select = !font_select;
      return;

   case '0':  case '1':  case '2':  case '3':  case '4':
   case '5':  case '6':  case '7':  case '8':
      if (!font_select) 
         ind_font = item - '0';
      else
         title_font = item - '0';
      break;

   default:
      break;
   }
   glutPostRedisplay();
   return;
}
/* End of menu */


/*
 *--------------------------------------------------------------------
 *
 *  motion:  Callback triggered by mouse motion with one or more mouse
 *           buttons pressed.  
 *
 *  If mode = CONSTRUCT_R, this function, operating in XOR write mode,
 *  draws the line segment with endpoints (xv1,yv1) and (xv2,yv2) 
 *  (thus erasing the previously drawn line segments), updates (xv2,
 *  yv2) to the object coordinates to the mouse position, and draws 
 *  the updated line segments.
 *
 *  If mode = ALTER_R, this function, operating in XOR write mode, 
 *  draws the two line segments with endpoints (xv1,yv1), (xv2,yv2),
 *  and (xv3,yv3) (thus erasing the previously drawn line segments), 
 *  updates (xv2,yv2) to the object coordinates to the mouse position, 
 *  and draws the updated line segments.
 *
 *  If mode = ALTER_VIEW_VOLUME, this function, operating in XOR write
 *  mode, draws an outline of the rectangle [xv1,xv2] X [yv1,yv2] 
 *  (thus erasing the previously drawn rectangle), updates (xv2,yv2)
 *  to the object coordinates corresponding to the mouse position, and
 *  draws an outline of the updated rectangle.
 *
 *  If mode = SELECT_CIRCLE, this function, operating in XOR write 
 *  mode, draws a circle of radius cr centered at (cx,cy), thus 
 *  erasing the previously drawn circle, changes cr to the distance
 *  from (cx,cy) to the mouse position, and draws the circle with the
 *  new radius.
 *
 *  If nnbv > 1 (and mode = ALTER_DOMAIN), this function, operating in
 *  XOR write mode, draws the set of line segments from (xv1,yv1) to 
 *  the neighbors of vertex kv1, erasing the previously drawn line 
 *  segments, updates (xv1,yv1) to the mouse position, and draws the 
 *  new set of line segments.  The indices of the neighbors of kv1 
 *  are stored in the first nnbv positions of lste.
 *
 *  glmesh2 functions required:  displayPosition, drawCircle, xlate
 *
 *--------------------------------------------------------------------
 */
static void motion(int x, int y)
{
   GLdouble xm, ym, zm;      /* Computed object coordinates */
   GLdouble x2, y2;
   GLint i, kv, nnb;

   xlate(x,y,0.0, &xm,&ym,&zm);
   if (mode == CONSTRUCT_R) {
/*
 *  Erase line segment from (xv1,yv1) to (xv2,yv2).
 */
      glBegin(GL_LINES);
        glVertex2d(xv1, yv1);
        glVertex2d(xv2, yv2);
      glEnd();
/*
 *  Update (xv2,yv2).
 */
      xv2 = xm;
      yv2 = ym;
/*
 *  Draw line segment from (xv1,yv1) to (xv2,yv2).
 */
      glBegin(GL_LINES);
        glVertex2d(xv1, yv1);
        glVertex2d(xv2, yv2);
      glEnd();
   } else if (mode == ALTER_R) {
/*
 *  Erase line segments from (xv1,yv1) to (xv2,yv2) and (xv2,yv2) to
 *  (xv3,yv3).
 */
      glBegin(GL_LINE_STRIP);
        glVertex2d(xv1, yv1);
        glVertex2d(xv2, yv2);
        glVertex2d(xv3, yv3);
      glEnd();
/*
 *  Update (xv2,yv2).
 */
      xv2 = xm;
      yv2 = ym;
/*
 *  Draw line segments from (xv1,yv1) to (xv2,yv2) and (xv2,yv2) to
 *  (xv3,yv3).
 */
      glBegin(GL_LINE_STRIP);
         glVertex2d(xv1, yv1);
         glVertex2d(xv2, yv2);
         glVertex2d(xv3, yv3);
      glEnd();
   } else if (mode == ALTER_VIEW_VOLUME) {
/*
 *  Erase rectangle [xv1,xv2] X [yv1,yv2].
 */
      glBegin(GL_LINE_LOOP);
        glVertex2d(xv1, yv1);
        glVertex2d(xv2, yv1);
        glVertex2d(xv2, yv2);
        glVertex2d(xv1, yv2);
      glEnd();
/*
 *  Update (xv2,yv2).
 */
      xv2 = xm;
      yv2 = ym;
/*
 *  Draw rectangle [xv1,xv2] X [yv1,yv2].
 */
      glBegin(GL_LINE_LOOP);
        glVertex2d(xv1, yv1);
        glVertex2d(xv2, yv1);
        glVertex2d(xv2, yv2);
        glVertex2d(xv1, yv2);
      glEnd();
   } else if (mode == SELECT_CIRCLE) {
/*
 *  mode = SELECT_CIRCLE.
 */
      drawCircle();
      cr = sqrt( (xm-cx)*(xm-cx) + (ym-cy)*(ym-cy) );
      drawCircle();
   } else if (nnbv > 1) {
/*
 *  mode = ALTER_DOMAIN and nnbv > 1.
 */
      nnb = nnbv;
      glBegin(GL_LINES);
      for (i = 0; i < nnb; i++) {
         kv = lste[i];
         x2 = vtxy[kv][0];
         y2 = vtxy[kv][1];
        glVertex2d(xv1, yv1);
        glVertex2d(x2, y2);
      }
      glEnd();
      xv1 = xm;
      yv1 = ym;
      glBegin(GL_LINES);
      for (i = 0; i < nnb; i++) {
         kv = lste[i];
         x2 = vtxy[kv][0];
         y2 = vtxy[kv][1];
        glVertex2d(xv1, yv1);
        glVertex2d(x2, y2);
      }
      glEnd();
   }
/*
 *  Display current mouse position in the upper left corner.
 */
   displayPosition(xm,ym);
   return;
}
/* End of motion */


/*
 *--------------------------------------------------------------------
 *
 *  mouse:  Callback triggered by a mouse button press or release with
 *          the mouse position in the primary window (except for 
 *          buttons attached to menus).
 *
 *  A left button press may be used to select and drag a vertex of the
 *  polygon boundary R, insert a new vertex in R by selecting a 
 *  boundary point other than a vertex, or delete a vertex by
 *  selecting it with the shift key engaged.
 *
 *  glmesh2 functions called:  intsec, plDistance, xlate
 *
 *--------------------------------------------------------------------
 */
static void mouse(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;      /* Computed object coordinates */
   GLdouble tols;
   GLint i, ip1, k, km1, km2, kp1;
   GLboolean intr;

   if (button != GLUT_LEFT_BUTTON  ||  np == 0) return;

   if (state == GLUT_DOWN) {
/*
 *  Compute the object coordinates (wx,wy) of the mouse cursor.
 */
      xlate(x,y,0.0, &wx,&wy,&wz);
/*
 *  Find a vertex R[k], if any, within squared distance tols of the
 *  mouse position (wx,wy).
 */
      tols = (sx <= sy) ? sx : sy;   /* tols = min(sx, sy)  */
      tols = 25.0*tols*tols;         /* 5-pixel tolerance   */
      intr = 0;
      for (k = 0; k <= np-1; k++) {
         if ((intr = ((wx-R[k][0])*(wx-R[k][0]) + 
                      (wy-R[k][1])*(wy-R[k][1]) <= tols))) break;
      }
      if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) {
         if (!intr  ||  np < 4) {
            return;
         } else {
/*
 *  Vertex selected with shift key down, and np > 3:  delete vertex k.
 */
         for (i = k; i < np-1; i++) {
            R[i][0] = R[i+1][0];
            R[i][1] = R[i+1][1];
         }
         np--;
         glutPostRedisplay();
         return;
         }
      }
      if (!intr) {
/*
 *  No shift key and no vertex selected.  Test for the mouse position
 *  within a tolerance of the polygon boundary.
 */
         tols = 0.64*tols;    /* 4-pixel tolerance  */
         intr = 0;
         for (k = 0; k < np; k++) {
            kp1 = (k+1)%np;
            if ((intr = plDistance(wx, wy, R[k][0], R[k][1], 
                                  R[kp1][0], R[kp1][1], tols))) break;
         }
         if (intr) {
/*
 *  Insert a new vertex after reallocating memory for R if necessary.
 */
            if (np >= npmax) {
               npmax *= 2;
               R = realloc(R, npmax*sizeof *R);
               if (R == NULL) {
                  printf("glmesh2:  Unable to allocate sufficient "
                         "memory for R.\n");
                  exit(1);
               }
            }
            for (i = np-1; i > k; i--) {
               R[i+1][0] = R[i][0];
               R[i+1][1] = R[i][1];
            }
            k++;
            R[k][0] = wx;
            R[k][1] = wy;
            np++;
         } else {
            return;
         }
      }
/*
 *  Set up for dragging (Function motion):  store R[k-1], R[k], and 
 *  R[k+1] in (xv1,yv1), (xv2,yv2), and (xv3,yv3), set XOR write mode,
 *  set newR = 1, mode = ALTER_R, set indxR = k, and disable the cursor.
 */
      km1 = (k-1+np)%np;
      xv1 = R[km1][0];
      yv1 = R[km1][1];
      xv2 = R[k][0];
      yv2 = R[k][1];
      kp1 = (k+1)%np;
      xv3 = R[kp1][0];
      yv3 = R[kp1][1];
      glEnable(GL_COLOR_LOGIC_OP);
      glLogicOp(GL_XOR);
      newR = 1;
      mode = ALTER_R;
      indxR = k;
      glutSetCursor(GLUT_CURSOR_NONE);
   } else {
/*
 *  Left button release:  test for an intersection.
 */
      if (mode != ALTER_R) return;
      k = indxR;
      km1 = (k-1+np)%np;
      km2 = (k-2+np)%np;
      intr = 0;
      for (i = 0; i < np; i++) {
         ip1 = (i+1)%np;
         if (i == km2  ||  i == km1  ||  i == k) continue;
         if (intsec(xv1,yv1,xv2,yv2,R[i][0],R[i][1],R[ip1][0],
                    R[ip1][1])) {
            intr = 1;
            break;
         }
      }
      if (!intr  &&  np > 3) {
         intr = intsec(xv2,yv2,xv3,yv3,R[km2][0],R[km2][1],R[km1][0],
                       R[km1][1]);
      }
      if (intr) {
/*
 *  Invalid choice:  erase the line segments.
 */
         glBegin(GL_LINE_STRIP);
            glVertex2d(xv1,yv1);
            glVertex2d(xv2,yv2);
            glVertex2d(xv3,yv3);
         glEnd();
         glutSwapBuffers();
      } else {
/*
 *  Update R.
 */
         R[k][0] = xv2;
         R[k][1] = yv2;
      }
/*
 *  Restore the inherited cursor, and disable XOR write mode.
 */
      mode = DEFAULT;
      glutSetCursor(GLUT_CURSOR_INHERIT);
      glDisable(GL_COLOR_LOGIC_OP);
      glutPostRedisplay();
   }
   return;
}
/* End of mouse */


/*
 *--------------------------------------------------------------------
 *
 *  moveRn:  Project all non-removable vertices in R (that can be
 *           moved without destroying the mesh) onto the current
 *           circle (if any).
 *
 *  glmesh2 functions called:  bbox, initR, intsec
 *
 *--------------------------------------------------------------------
 */
static void moveRn(void)
{
   double s, x1, x2, x3, xp, y1, y2, y3, yp;
   int i, k, kt, kv, kv1, kv2, kv3 = 0, kvn;
   unsigned char intr, visible;
/*
 *  Test for a circle.
 */
   if (cr <= 0.0) return;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
/*
 *  Outer loop on non-removable vertices kv1 = lnrv[k] in R.
 */
   for (k = 0; k < nn; k++) {
      kv1 = lnrv[k];
      if (!listf[kv1]) continue;
/*
 *  Compute the coordinates (xp,yp) of the projection p of kv1
 *  onto the circle.  If kv1 is not at the center (cx,cy), it can
 *  be projeced onto the circle (but it could end up anywhere on
 *  the circle if it starts too near the center).
 */
      x1 = vtxy[kv1][0];
      y1 = vtxy[kv1][1];
      s = sqrt( (x1-cx)*(x1-cx) + (y1-cy)*(y1-cy) );
      if (s <= 0.0) continue;
      s = cr/s;
      xp = cx + s*(x1-cx);
      yp = cy + s*(y1-cy);
/*
 *  The new point p is visible iff p Left kv2 -> kv3 for all pairs of
 *  adjacent neighbors (kv2,kv3) of kv1, and there is no intersection
 *  of the path p-kv1 with a boundary edge not containing kv1 in the
 *  same boundary curve.
 */
      visible = 1;
      kt = lvtx[kv1];
      while (kt >= 0) {
         if (ltri[kt][0] == kv1)
            i = 1;
         else if (ltri[kt][1] == kv1)
            i = 2;
         else
            i = 0;
         kv2 = ltri[kt][i];
         kv3 = ltri[kt][(i+1)%3];
         x2 = vtxy[kv2][0];
         y2 = vtxy[kv2][1];
         x3 = vtxy[kv3][0];
         y3 = vtxy[kv3][1];
         if ( (x3-x2)*(yp-y2) < (y3-y2)*(xp-x2) ) {
            visible = 0;
            break;   
         }
         kt = ltri[kt][i+3];
      }
      if (visible) {
/*
 *  Set kv2 and kv3 to the first and last neighbors of kv1, and 
 *  test for intersections.
 */
         intr = 0;
         kt = lvtx[kv1];
         if (ltri[kt][0] == kv1)
            i = 1;
         else if (ltri[kt][1] == kv1)
            i = 2;
         else
            i = 0;
         kv2 = ltri[kt][i];
         kv = kv2;
         while (kv != kv3) {
            kt = lvtx[kv];
            if (ltri[kt][0] == kv)
               i = 1;
            else if (ltri[kt][1] == kv)
               i = 2;
            else
               i = 0;
            kvn = ltri[kt][i];
            if (intsec(xp,yp,x1,y1,vtxy[kv][0],vtxy[kv][1],
                       vtxy[kvn][0],vtxy[kvn][1])) {
               intr = 1;
               break;
            }
            kv = kvn;
         }
         if (intr) visible = 0;
      }
      if (visible) {
/*
 *  Move kv1 to its new position.
 */
         vtxy[kv1][0] = xp;
         vtxy[kv1][1] = yp;
      }
   }
/*
 *  Update the bounding box.
 */
   bbox();
   return;
}
/* End of moveRn */


/*
 *--------------------------------------------------------------------
 *
 *  optimizeRt:  Optimize R with Delaunay edge swaps in convex    
 *               quadrilaterals consisting of pairs of adjacent 
 *               triangles in R.
 *
 *  The indices of the interior edges associated with pairs of 
 *  adjacent triangles contained in R are stored in lste.  Then a
 *  sequence of sweeps is executed in which each sweep consists of
 *  applying the swap test (empty circumcircle test) and appropriate
 *  swaps to the interior edges in the order in which they are stored.
 *  The iteration is repeated until no swap occurs or maxit iterations
 *  have been performed.  The bound on the number of iterations 
 *  (defined below) may be necessary to prevent an infinite loop 
 *  caused by cycling (reversing the effect of a previous swap) due to
 *  floating point inaccuracy when four or more vertices are nearly 
 *  cocircular.
 *
 *  glmesh2 functions called:  initR, swap, swptst
 *
 *--------------------------------------------------------------------
 */
static void optimizeRt(void)
{
   int maxit = 4;          /* Maximum number of iterations */
   unsigned char swp;      /* Flag set by a swap */
   int in1, in2, io1, io2;
   int i, iter, j, ke, kt, ktn, le;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
/*
 *  Store the interior edge indices in lste with count le.
 *  Each pair of adjacent triangles in R is encountered twice
 *  (once for each triangle).  The corresponding edge is only
 *  stored once by associating it with the smaller triangle index.
 */
   le = 0;
   for (j = 0; j < ntr; j++) {
      kt = listtr[j];
      for (i = 0; i < 3; i++) {
         ktn = ltri[kt][i+3];
         if (ktn < kt) continue;
         if (listf[ltri[ktn][0]]  &&  listf[ltri[ktn][1]]  &&
             listf[ltri[ktn][2]]) {
            lste[le] = ltri[kt][i+6];
            le++;
         }
      }
   }
   if (le == 0) return;
/*
 *  Top of outer loop.
 */
   for (iter = 0; iter < maxit; iter++) {
      swp = 0;
/*
 *  Inner loop on interior edges.
 */
      for (j = 0; j < le; j++) {
         ke = lste[j];
         kt = ledg[ke];
         if (ltri[kt][6] == ke)
            i = 0;
         else if (ltri[kt][7] == ke)
            i = 1;
         else
            i = 2;
         in1 = ltri[kt][i];
         io1 = ltri[kt][(i+1)%3];
         io2 = ltri[kt][(i+2)%3];
         kt = ltri[kt][i+3];
         if (ltri[kt][6] == ke)
            i = 0;
         else if (ltri[kt][7] == ke)
            i = 1;
         else
            i = 2;
         in2 = ltri[kt][i];
         if (swptst(in1,in2,io1,io2)) {
            swap(ke, 0);
            swp++;
         }
      }
      if (!swp) break;
   }
   return;
}
/* End of optimizeRt */


/*
 *--------------------------------------------------------------------
 *
 *  optimizeRv:  Optimize R by applying a single iteration of a CVT
 *               (Centroid Voronoi Tessellation) smoothing operation:
 *               sweep through the removable vertices, moving each 
 *               boundary vertex to the midpoint of its predecessor
 *               and successor, and each interior vertex to a weighted
 *               sum (convex combination) of the centroids of the 
 *               triangles that contain the vertex, where the weights
 *               are the triangle areas relative to the patch area.  
 *               This is similar to Laplacian smoothing by a Gauss-
 *               Seidel iteration.  If global variable bdry_layer =
 *               TRUE, each interior vertex adjacent to a boundary 
 *               vertex is moved to a location that depends only on
 *               the triangles that contain the vertex and contain a
 *               boundary vertex, so that the new location is usually
 *               closer to the boundary.
 *
 *  glmesh2 function called:  initR
 *
 *--------------------------------------------------------------------
 */
static void optimizeRv(void)
{
   double s, sumw, w, x, x0, x1, x2, y, y0, y1, y2;
   unsigned char bdry, bdry1, bdry2, removable, visible;
   int i, j, k, kt, kt1, kv0, kv1, kv2, kvn, nnb;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
/*
 *  Loop on vertices kv0 in R.
 */
   for (k = 0; k < nvr; k++) {
      kv0 = listvr[k];
      x0 = vtxy[kv0][0];
      y0 = vtxy[kv0][1];
/*
 *  Test for kv0 removable.
 */
      removable = 1;
      for (i = 0; i < nn; i++) {
         if (lnrv[i] == kv0) {
            removable = 0;
            break;
         }
      }
      if (!removable) continue;
/*
 *  Store the neighbors of kv0 in the first nnb positions of lste,
 *  where the first neighbor is also the last neighbor (stored twice)
 *  unless kv0 is a boundary vertex.
 */
      nnb = 0;
      kt1 = lvtx[kv0];
      kt = kt1;
      do {
         if (ltri[kt][0] == kv0)
            i = 1;
         else if (ltri[kt][1] == kv0)
            i = 2;
         else
            i = 0;
         lste[nnb] = ltri[kt][i];
         nnb++;
         kvn = ltri[kt][(i+1)%3];
         kt = ltri[kt][i+3];
      } while (kt >= 0  &&  kt != kt1);
      lste[nnb] = kvn;
      nnb++;
      if (kt < 0) {
/*
 *  Boundary vertex:  compute the new position (x,y) as the midpoint 
 *  of the first and last neighbors kv1 and kv2.
 */
         kv1 = lste[0];
         kv2 = lste[nnb-1];
         x = 0.5*(vtxy[kv1][0] + vtxy[kv2][0]);
         y = 0.5*(vtxy[kv1][1] + vtxy[kv2][1]);
         visible = 1;
      } else {
/*
 *  Interior vertex:  set bdry = TRUE iff bdry_layer = TRUE and a 
 *  neighbor kv1 of kv0 is a boundary vertex.
 */
	 bdry = 0;
         if (bdry_layer) {
            for (i = 0; i < nnb-1; i++) {
               kv1 = lste[i];
               kt = lvtx[kv1];
               if (ltri[kt][0] == kv1)
                  j = 2;
               else if (ltri[kt][1] == kv1)
                  j = 0;
               else
                  j = 1;
               if ((bdry = (ltri[kt][j+3] < 0))) break; 
            }
         }
/*
 *  Compute the potential new position as a weigted sum (x,y) of the
 *  centroids of the triangles (kv0,kv1,kv2).
 *
 *  Initialization:
 */
         x = 0.0;
         y = 0.0;
         sumw = 0.0;
	 bdry1 = 0;
	 bdry2 = 0;
         if (bdry) {
            kv1 = lste[0];
            kt = lvtx[kv1];
            if (ltri[kt][0] == kv1)
               j = 2;
            else if (ltri[kt][1] == kv1)
               j = 0;
            else
               j = 1;
            bdry2 = (ltri[kt][j+3] < 0); 
	 }
         for (i = 0; i < nnb-1; i++) {
            bdry1 = bdry2;
            kv1 = lste[i];
            kv2 = lste[i+1];
            if (bdry) {
/*
 *  Bypass the triangle if neither kv1 nor kv2 is a boundary vertex.
 */
               kt = lvtx[kv2];
               if (ltri[kt][0] == kv2)
                  j = 2;
               else if (ltri[kt][1] == kv2)
                  j = 0;
               else
                  j = 1;
               bdry2 = (ltri[kt][j+3] < 0); 
	       if (!bdry1  &&  !bdry2) continue;
            }
            x1 = vtxy[kv1][0];
            y1 = vtxy[kv1][1];
            x2 = vtxy[kv2][0];
            y2 = vtxy[kv2][1];
            w = (x2-x1)*(y0-y1) - (y2-y1)*(x0-x1);
            sumw += w;
            x += w*(x0 + x1 + x2);
            y += w*(y0 + y1 + y2);
         }
         s = 1.0/(3.0*sumw);
         x *= s;
         y *= s;
/*
 *  Test for (x,y) visible from the neighbors:  (x,y) Left kv1 -> kv2
 *  for all pairs of adjacent neighbors (kv1,kv2).
 */
         visible = 1;
         for (i = 0; i < nnb-1; i++) {
            kv1 = lste[i];
            kv2 = lste[i+1];
            x1 = vtxy[kv1][0];
            y1 = vtxy[kv1][1];
            x2 = vtxy[kv2][0];
            y2 = vtxy[kv2][1];
            if ( (x2-x1)*(y-y1) < (y2-y1)*(x-x1) ) {
               visible = 0;
               break;   
            }
         }
      }
      if (!visible) continue;
/*
 *  Move kv0 to its new position.
 */
      vtxy[kv0][0] = x;
      vtxy[kv0][1] = y;
   }
   return;
}
/* End of optimizeRv */


/*
 *--------------------------------------------------------------------
 *
 *  passiveMotion:  Callback triggered by mouse motion with no mouse
 *                  buttons pressed.  
 *
 *  This function displays the object coordinates of the mouse
 *  position.
 *
 *  glmesh2 functions called:  displayPosition, xlate
 *
 *--------------------------------------------------------------------
 */
static void passiveMotion(int x, int y)
{
   GLdouble xm, ym, zm;           /* Computed object coordinates */
/*
 *  Compute the object coordinates (xm,ym) of the mouse cursor.
 */
   xlate(x,y,0.0, &xm,&ym,&zm);
/*
 *  Display current mouse position in the upper left corner.
 */
   displayPosition(xm,ym);
   return;
}
/* End of passiveMotion */


/*
 *--------------------------------------------------------------------
 *
 *  plDistance:  Compute a point-line distance.  The return value is
 *               TRUE iff the projection p* of point p onto the line
 *               defined by p1 and p2 lies on the line segment p1-p2
 *               and satisfied <p*-p,p*-p> <= tols.
 *
 *  On input:  p = (xp,yp), p1 = (x1,y1), p2 = (x2,y2), and tols is
 *             the maximum squared distance between p and p* for which
 *             the return value is TRUE.
 *
 *  glmesh2 functions called:  None
 *  
 *--------------------------------------------------------------------
 */
static unsigned char plDistance(double xp, double yp, double x1, 
                                double y1, double x2, double y2, 
                                double tols)
{
   double cp, dp, ds12, dx12, dx1p, dy12, dy1p;
/*
 *  p* = p1 + s*(p2-p1), where <p-p*,p2-p1> = 0, so that
 *  s = dp/ds12 for dp = <p-p1,p2-p1> and ds12 = <p2-p1,p2-p1>.
 *  Then p* lies on the line segment iff 0 <= s <= 1 for ds12 > 0.
 *  The return value is FALSE if p1 = p2.
 */
   dx12 = x2-x1;
   dy12 = y2-y1;
   ds12 = dx12*dx12 + dy12*dy12;
   dx1p = xp-x1;
   dy1p = yp-y1;
   dp = dx1p*dx12 + dy1p*dy12;
   if (dp < 0.0  ||  dp > ds12) return 0;
/*
 *  The squared distance is ds = <p*-p,p*-p> = cp*cp/ds12 for 
 *  cross product (z-component) cp = (p2-p1) X (p-p1).
 */
   cp = dx12*dy1p - dy12*dx1p;
   return (cp*cp <= tols*ds12);
}
/* End of plDistance */


/*
 *--------------------------------------------------------------------
 *
 *  refineRe:  Refine R by adding a new vertex at the midpoint of each
 *             edge in R.  The new vertex is connected to the end-
 *             points of the original edge and to the one or two 
 *             vertices opposite the original edge, thus replacing 
 *             each triangle involved by two new triangles.  Each 
 *             triangle contained in R is partitioned into four sub-
 *             triangles by three interior edges.  Following the 
 *             additions, an edge swap is applied to the first 
 *             interior edge that was added to each triangle in R so
 *             that, if the original triangle was equilateral, the
 *             four new triangles are also equilateral.
 *
 *  glmesh2 functions called:  initR, split, swap
 *
 *--------------------------------------------------------------------
 */
static void refineRe(void)
{
   unsigned char ktinR = 0, ktninR = 0;
   int i, k, ke, ken, kt, ktn, kv0, kv3, le, ne0, ner0, nv0;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
   if (!ner) return;
/*
 *  Push nv onto rstack.
 */
   if (nref == nrefmax) {
      rstack = realloc (rstack, nrefmax*sizeof *rstack);
      if (rstack == NULL) {
         printf("glmesh2:  Unable to allocate sufficient memory "
                "for rstack.\n");
         exit(1);
      }
   }
   rstack[nref] = nv;
   nref++;
   if (nv + ner > nvmax) {
/*
 *  Increase nvmax to twice the new nv value, and re-allocate storage 
 *  for vtxy, uv, lvtx, ltri, ledg, lnrv, lcrv, listf, lister, listtr, 
 *  listvr, and lste.
 */
      nvmax = 2*(nv+ner);
      vtxy = realloc (vtxy, nvmax*sizeof *vtxy);
      uv = realloc (uv, nvmax*sizeof *uv);
      lvtx = realloc (lvtx, nvmax*sizeof *lvtx);
      ltri = realloc (ltri, 2*nvmax*sizeof *ltri);
      ledg = realloc (ledg, 3*nvmax*sizeof *ledg);
      lnrv = realloc (lnrv, nvmax*sizeof *lnrv);
      lcrv = realloc (lcrv, (nvmax/3)*sizeof *lcrv);
      listf = realloc (listf, nvmax*sizeof *listf);
      lister = realloc (lister, 3*nvmax*sizeof *lister);
      listtr = realloc (listtr, 2*nvmax*sizeof *listtr);
      listvr = realloc (listvr, nvmax*sizeof *listvr);
      lste = realloc (lste, 3*nvmax*sizeof *lste);
      if (vtxy == NULL  ||  uv == NULL  ||  lvtx == NULL  ||  
          ltri == NULL  ||  ledg == NULL  ||  lnrv == NULL  ||  
          lcrv == NULL  ||  listf == NULL  ||  lister == NULL  ||  
          listtr == NULL  ||  listvr == NULL  ||  lste == NULL) {
         printf("glmesh2:  Unable to allocate sufficient memory.\n");
         exit(1);
      }
   }
/*
 *  Initialization:  
 *
 *  le = Number of elements in lste,
 *  ner0 = Initial length of lister,
 *  ne0 = Initial edge count,
 *  nv0 = Initial vertex count.
 */
   le = 0;
   ner0 = ner;
   ne0 = ne;
   nv0 = nv;
/*
 *  Loop on edges ke in R.
 */
   for (k = 0; k < ner0; k++) {
      ke = lister[k];
      split(ke);
/*
 *  The first edges added to the interior of the triangles that
 *  were split have endpoint indices nv-1 and kv0 < nv0.  If either
 *  triangle is in R, the corresponding edge is added to the list
 *  of edges lste to be swapped.
 */
      kt = ledg[ke];
      if (ltri[kt][6] == ke)
         i = 0;
      else if (ltri[kt][7] == ke)
         i = 1;
      else
         i = 2;
      kv0 = ltri[kt][i];
      ktn = ltri[kt][i+3];
      ktinR = listf[kv0];
      if (ktinR  &&  kv0 < nv0) {
         ken = ne-1;
         if (ktn >= 0) ken--;
         lste[le] = ken;
         le++;
      }
/*
 *  The first edge added to the interior of triangle ktn, if ktn is
 *  in R, must also be added to the list.
 */
      if (ktn >= 0) {
         if (ltri[ktn][3] == kt)
            i = 0;
         else if (ltri[ktn][4] == kt)
            i = 1;
         else
            i = 2;
         kv3 = ltri[ktn][i];
         ktninR = listf[kv3];
         if (ktninR  &&  kv3 < nv0) {
            lste[le] = ne-1;
            le++;
         }
      }
/*
 *  Update lists of elements in R:  one new vertex, zero to two new
 *  triangles, and one to three new edges.
 */
      listvr[nvr] = nv-1;
      nvr++;
      listf[nv-1] = 1;
      if (ktn >= 0) {
         lister[ner] = ne-3;
         ner++;
         if (ktinR) {
            listtr[ntr] = nt-2;
            ntr++;
            lister[ner] = ne-2;
            ner++;
         }
         if (ktninR) {
            listtr[ntr] = nt-1;
            ntr++;
            lister[ner] = ne-1;
            ner++;
         }
      } else {
         lister[ner] = ne-2;
         ner++;
         if (ktinR) {
            listtr[ntr] = nt-1;
            ntr++;
            lister[ner] = ne-1;
            ner++;
         }
      }
   }  /* end of loop on edges */
/*
 *  Apply swaps to the edges in lste.
 */
   for (k = 0; k < le; k++) {
      ke = lste[k];
      swap(ke, 0);
   }
   return;
}
/* End of refineRe */


/*
 *--------------------------------------------------------------------
 *
 *  refineRt:  Refine R by adding triangle barycenters.  Refer to
 *             Function insertb.
 *
 *  glmesh2 functions called:  initR, insertb
 *
 *--------------------------------------------------------------------
 */
static void refineRt(void)
{
   int k, ntr0;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
   if (!ntr) return;
/*
 *  Push nv onto rstack.
 */
   if (nref == nrefmax) {
      rstack = realloc (rstack, nrefmax*sizeof *rstack);
      if (rstack == NULL) {
         printf("glmesh2:  Unable to allocate sufficient memory "
                "for rstack.\n");
         exit(1);
      }
   }
   rstack[nref] = nv;
   nref++;
   if (nv + ntr > nvmax) {
/*
 *  Increase nvmax to twice the new nv value, and re-allocate storage 
 *  for vtxy, uv, lvtx, ltri, ledg, lnrv, lcrv, listf, lister, listtr, 
 *  listvr, and lste.
 */
      nvmax = 2*(nv+ntr);
      vtxy = realloc (vtxy, nvmax*sizeof *vtxy);
      uv = realloc (uv, nvmax*sizeof *uv);
      lvtx = realloc (lvtx, nvmax*sizeof *lvtx);
      ltri = realloc (ltri, 2*nvmax*sizeof *ltri);
      ledg = realloc (ledg, 3*nvmax*sizeof *ledg);
      lnrv = realloc (lnrv, nvmax*sizeof *lnrv);
      lcrv = realloc (lcrv, (nvmax/3)*sizeof *lcrv);
      listf = realloc (listf, nvmax*sizeof *listf);
      lister = realloc (lister, 3*nvmax*sizeof *lister);
      listtr = realloc (listtr, 2*nvmax*sizeof *listtr);
      listvr = realloc (listvr, nvmax*sizeof *listvr);
      lste = realloc (lste, 3*nvmax*sizeof *lste);
      if (vtxy == NULL  ||  uv == NULL  ||  lvtx == NULL  ||  
          ltri == NULL  ||  ledg == NULL  ||  lnrv == NULL  ||  
          lcrv == NULL  ||  listf == NULL  ||  lister == NULL  ||  
          listtr == NULL  ||  listvr == NULL  ||  lste == NULL) {
         printf("glmesh2:  Unable to allocate sufficient memory.\n");
         exit(1);
      }
   }
/*
 *  Loop on triangles in R, inserting the barycenters.
 */
   ntr0 = ntr;
   for (k = 0; k < ntr0; k++) {
      insertb(listtr[k]);
   }
   return;
}
/* End of refineRt */


/*
 *--------------------------------------------------------------------
 *
 *  removeRt:  Remove the triangles in R from the triangle mesh (and
 *             from R), creating new boundary curves if necessary.  It
 *             may not be possible to remove all triangles if removal
 *             would render the mesh invalid.
 *
 *  glmesh2 functions called:  deltri, initR
 *
 *--------------------------------------------------------------------
 */
static void removeRt(void)
{
   int btest = 1;    /* Flag specifying a test for a boundary edge */
   int k;            /* Index for listtr */
   int kt0;          /* Triangle index */
   int nd = 1;       /* Flag specifying another sweep required */
   int nv0 = nv;     /* Initial vertex count */
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
   
   while (nd) {
      nd = 0;
      k = 0;
      while (k < ntr) {
         kt0 = listtr[k];
         if (!btest  ||  (ltri[kt0][3] < 0  ||
             ltri[kt0][4] < 0  ||  ltri[kt0][5] < 0)) {
            if (deltri(kt0)) {
               nd++;
               if (!btest) btest = 1;
            } else {
               k++;
            } 
         } else {
            k++;
         }
      }
      if (!nd  &&  btest) {
         btest = 0;
         nd = 1;
      }
   }
/*
 *  The refinement stack is emptied if a vertex was deleted.
 */
   if (nv < nv0) nref = 0;
   return;
}
/* End of removeRt */


/*
 *--------------------------------------------------------------------
 *
 *  reshape:  Set up a viewport, construct a projection, and 
 *            initialize the modelview matrix.
 *
 *  Global variables required:  aspect, dx, dy, sx, sy, xc, yc
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void reshape(int w, int h)
{
   GLdouble r;                   /* Aspect ratio of the view volume */
   int imax, imin, jmax, jmin;   /* Window coordinates of viewport */
   int ni, nj;                   /* Width and height of viewport */
   if (aspect) {
/*
 *  Compute viewport corner coordinates such that the aspect ratio r
 *  of the view volume is preserved (assuming square pixels) and the
 *  viewport is centered in the window and as large as possible.
 */
      r = dx/dy;
      if (r <= (GLdouble) w/(GLdouble) h) {
         imin = (int) ((w-r*h)/2.0);
         imax = imin + (int) (r*h+0.5);
         if (imax > w) imax = w;
         jmin = 0;
         jmax = h;
      } else {
         imin = 0;
         imax = w;
         jmin = (int) ((h-w/r)/2.0);
         jmax = jmin + (int) (w/r+0.5);
      }
   } else {
      imin = 0;
      imax = w;
      jmin = 0;
      jmax = h;
   }
   ni = imax-imin;
   nj = jmax-jmin;
/*
 *  Store global variable sx and sy, and create the viewport mapping.
 */
   sx = dx/(GLdouble) ni;
   sy = dy/(GLdouble) nj;

   glViewport(imin, jmin, ni, nj);

   glMatrixMode(GL_PROJECTION);
   glLoadIdentity();
/*
 *  Orthographic projection.
 */
   gluOrtho2D(xc-dx/2.0, xc+dx/2.0, yc-dy/2.0, yc+dy/2.0);

   glMatrixMode(GL_MODELVIEW);
   glLoadIdentity();

   return;
}
/* End of reshape */


/*
 *--------------------------------------------------------------------
 *
 *  specialKey:  Move center (xc,yc) in response to arrow keys, zoom
 *               in or out on Page Up or Page Down keys, or restore
 *               default view volume if Home key is pressed.
 *
 *  Global variables required:  dx, dy, xc, yc, xmin, ymin, xmax, ymax
 *
 *  glmesh2 function called:  reshape
 *
 *--------------------------------------------------------------------
 */
static void specialKey(int key, int x, int y)
{
   switch (key) {
      case GLUT_KEY_UP:
         if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) {
            yc += 0.25*(ymax-ymin);
         } else {
            yc += 0.05*(ymax-ymin);
         }
         break;
      case GLUT_KEY_DOWN:
         if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) {
            yc -= 0.25*(ymax-ymin);
         } else {
            yc -= 0.05*(ymax-ymin);
         }
         break;
      case GLUT_KEY_RIGHT:
         if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) {
            xc += 0.25*(xmax-xmin);
         } else {
            xc += 0.05*(xmax-xmin);
         }
         break;
      case GLUT_KEY_LEFT:
         if (glutGetModifiers() == GLUT_ACTIVE_SHIFT) {
            xc -= 0.25*(xmax-xmin);
         } else {
            xc -= 0.05*(xmax-xmin);
         }
         break;

      case GLUT_KEY_PAGE_UP:     /* Zoom in */
         dx *= 0.8;
         dy *= 0.8;
         break;
      case GLUT_KEY_PAGE_DOWN:   /* Zoom out */
         dx *= 1.25;
         dy *= 1.25;
         break;

      case GLUT_KEY_HOME:        /* Restore default view volume */
         dx = 1.2*(xmax-xmin);
         dy = 1.2*(ymax-ymin);
         xc = (xmin+xmax)/2.0;
         yc = (ymin+ymax)/2.0;
         break;

      default:
         break;
   }
   reshape(glutGet(GLUT_WINDOW_WIDTH),
           glutGet(GLUT_WINDOW_HEIGHT));
   glutPostRedisplay();

   return;
}
/* End of specialKey */


/*
 *--------------------------------------------------------------------
 *
 *  split:  Add a new vertex at the midpoint of an edge, splitting the
 *          edge into a pair of edges, and replacing each triangle 
 *          containing the edge by two triangles.  The uv values at
 *          the midpoint are obtained by linear interpolation.
 *
 *  On input:  ke = index of the edge to be split.
 *
 *  It is necessary that nv < nvmax so that there is sufficient room
 *  in arrays vtxy, ltri, lvtx, and ledg to add another vertex, two
 *  more triangles, and three more edges.  No variables are altered
 *  if nv >= nvmax.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void split(int ke)
{
   int i, j, ke1, ke2 = 0, kt, kt1, kt2, ktn, kv0, kv1, kv2, kv3 = 0;
   if (nv >= nvmax) return;
/*
 *  Set (kv0,kv1,kv2) to the vertex indices of the triangle kt that
 *  contains ke, where kv0 = ltri[kt0][i] and ke = ltri[kt0][i+6].
 */
   kt = ledg[ke];
   if (ltri[kt][6] == ke)
      i = 0;
   else if (ltri[kt][7] == ke)
      i = 1;
   else
      i = 2;
   kv0 = ltri[kt][i];
   kv1 = ltri[kt][(i+1)%3];
   kv2 = ltri[kt][(i+2)%3];
/*
 *  Append the new vertex.
 */
   vtxy[nv][0] = 0.5*(vtxy[kv1][0]+vtxy[kv2][0]);
   vtxy[nv][1] = 0.5*(vtxy[kv1][1]+vtxy[kv2][1]);
   uv[nv][0] = 0.5*(uv[kv1][0]+uv[kv2][0]);
   uv[nv][1] = 0.5*(uv[kv1][1]+uv[kv2][1]);
/*
 *  Set ktn to the index of the triangle opposite kv0, and set kt1
 *  and ke1 to the indices of the triangle and edge opposite kv1.
 */
   ktn = ltri[kt][i+3];
   kt1 = ltri[kt][(i+1)%3+3];
   ke1 = ltri[kt][(i+1)%3+6];
/*
 *  If i = j = 0, the initial ltri rows are
 *
 *    kt:   (kv0, kv1, kv2, ktn, kt1, x, ke, ke1, x)
 *    ktn:  (kv3, kv2, kv1, kt, x, kt2, ke, x, ke2) if ktn > 0,
 *
 *  and the new columns, following the split are
 *
 *    kt:    (kv0, kv1, nv, ktn, nt, x, ke, ne+1, x)
 *    nt:    (kv0, nv, kv2, nt+1, kt1, kt, ne, ke1, ne+1)
 *    nt+1:  (kv3, kv2, nv, nt, ktn, kt2, ne, ne+2, ke2)
 *    ktn:   (kv3, nv, kv1, kt, x, nt+1, ke, x, ne+2)
 * 
 *  where the entries labeled x retain their values and need not be
 *  referenced, nt+1 = -1 if ktn = -1, and the last two rows only 
 *  exist if ktn >= 0.  The changes implicit in the above lists are 
 *  identical for other cyclical orderings (other values of i and j).
 */
   ltri[kt][(i+2)%3] = nv;
   ltri[kt][(i+1)%3+3] = nt;
   ltri[kt][(i+1)%3+6] = ne+1;
   ltri[nt][0] = kv0;
   ltri[nt][1] = nv;
   ltri[nt][2] = kv2;
   ltri[nt][3] = -1;
   ltri[nt][4] = kt1;
   ltri[nt][5] = kt;
   ltri[nt][6] = ne;
   ltri[nt][7] = ke1;
   ltri[nt][8] = ne+1;
   kt2 = -1;
   if (ktn >= 0) {
/*
 *  Find j such such that kv3 = ltri[ktn][j] is the vertex opposite
 *  edge ke, and set kt2 and ke2 to the triangle and edge opposite 
 *  kv1 in ktn.
 */
      if (ltri[ktn][0] == kv1)
         j = 1;
      else if (ltri[ktn][1] == kv1)
         j = 2;
      else
         j = 0;
      kv3 = ltri[ktn][j];
      kt2 = ltri[ktn][(j+2)%3+3];
      ke2 = ltri[ktn][(j+2)%3+6];

      ltri[nt][3] = nt+1;
      ltri[ktn][(j+1)%3] = nv;
      ltri[ktn][(j+2)%3+3] = nt+1;
      ltri[ktn][(j+2)%3+6] = ne+2;

      ltri[nt+1][0] = kv3;
      ltri[nt+1][1] = kv2;
      ltri[nt+1][2] = nv;
      ltri[nt+1][3] = nt;
      ltri[nt+1][4] = ktn;
      ltri[nt+1][5] = kt2;
      ltri[nt+1][6] = ne;
      ltri[nt+1][7] = ne+2;
      ltri[nt+1][8] = ke2;
   }
/*
 *  Correct neighboring triangle ltri entries.
 */
   if (kt1 >= 0) {
      if (ltri[kt1][3] == kt)
         ltri[kt1][3] = nt;
      else if (ltri[kt1][4] == kt)
         ltri[kt1][4] = nt;
      else
         ltri[kt1][5] = nt;
   }
   if (kt2 >= 0) {
      if (ltri[kt2][3] == ktn)
         ltri[kt2][3] = nt+1;
      else if (ltri[kt2][4] == ktn)
         ltri[kt2][4] = nt+1;
      else
         ltri[kt2][5] = nt+1;
   }
/*  
 *  Update lvtx.
 */
   if (lvtx[kv2] == kt) lvtx[kv2] = nt;
   if (lvtx[kv2] == ktn) lvtx[kv2] = nt+1;
   if (ktn >= 0) {
      if (lvtx[kv3] == ktn) lvtx[kv3] = nt+1;
   }
   lvtx[nv] = nt;
/*
 *  Update ledg, preserving the property that the larger of the two 
 *  triangle indices is used.
 */
   ledg[ne] = nt;
   ledg[ne+1] = nt;
   ledg[ke1] = nt;
   if (ktn >= 0) {
      ledg[ne] = nt+1;
      ledg[ne+2] = nt+1;
      ledg[ke2] = nt+1;
   }
/*
 *  Update the counts.
 */
   nv++;
   if (ktn < 0) {
      nt++;
      ne += 2;
      nb++;
   } else {
      nt += 2;
      ne += 3;
   }
   return;
}
/* End of split */


/*
 *--------------------------------------------------------------------
 *
 *  storeV:  Compute and store global parameter values that depend on 
 *           the input data:  dx, dy, xc, xmax, xmin, yc, ymax, and
 *           ymin, and allocate storage for R and rstack.
 *
 *  Global variables required:
 *
 *      npmax = Initial size for R.
 *
 *      nrefmax = Initial size for rstack.
 *
 *      nv = Number of vertices in the triangle mesh.
 *
 *      R = Array dimensioned npmax by 2 containing the Cartesian 
 *          coordinates of a sequence of vertices defining the 
 *          boundary of the polygonal region.
 *
 *      rstack = Array of length nrefmax containing a sequence of
 *               vertex counts nv associated with mesh refinements.
 *
 *      vtxy = Array dimensioned nv by 2 containing the Cartesian 
 *             coordinates of the triangle mesh vertices.
 *
 *  glmesh2 function called:  bbox
 *
 *-------------------------------------------------------------------
 */
static void storeV(void)
{
/*
 *  Allocate storage.
 */
   R = malloc (npmax*sizeof *R);
   rstack = malloc (nrefmax*sizeof *rstack);
/*
 *  Compute the range of vertex coordinates (bounding box corner
 *  coordinates):  [xmin,xmax] X [ymin,ymax].
 */
   bbox();
/*
 *  Store the view volume dimensions and center.
 */
   dx = 1.2*(xmax-xmin);
   dy = 1.2*(ymax-ymin);
   xc = (xmin+xmax)/2.0;
   yc = (ymin+ymax)/2.0;
   return;
}
/* End of storeV */


/*
 *--------------------------------------------------------------------
 *
 *  swap:  Modify a triangle mesh by swapping diagonal edges in a con-
 *         vex quadrilateral defined by a pair of adjacent triangles 
 *         (triangles that share a common edge).  For a given set of 
 *         vertices, any triangle mesh can be transformed into any 
 *         other triangle mesh by a sequence of such swaps.  
 *
 *  On input:  
 *
 *      ke = Index of the edge to be swapped.
 *
 *      updateR = Flag with value 1 iff lister and listtr are to be
 *                updated.  This is not necessary if it is known
 *                that the triangles containing ke are in R.
 *
 *  It is assumed that the quadrilateral is convex.  No variables are
 *  altered if edge ke is not shared by two triangles.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void swap(int ke, int updateR)
{
   int i, j, k, ken1, ken2, kt1, kt2, ktn1, ktn2, kv1, kv2, kvn1, 
       kvn2;
/*
 *  Set kt1 and kt2 to the adjacent triangles with shared edge ke,
 *  set (kv1,kvn1,kvn2) and (kv2,kvn2,kvn1) to the vertex indices of 
 *  triangles kt1 and kt2, respectively, and set i and j to the column
 *  indices of kv1 in kt1 and kv2 in kt2, respectively.
 */
   kt1 = ledg[ke];
   if (ltri[kt1][6] == ke)
      i = 0;
   else if (ltri[kt1][7] == ke)
      i = 1;
   else
      i = 2;
   kt2 = ltri[kt1][i+3];
   if (kt2 < 0) return;
   if (ltri[kt2][6] == ke)
      j = 0;
   else if (ltri[kt2][7] == ke)
      j = 1;
   else
      j = 2;
   kv1 = ltri[kt1][i];
   kv2 = ltri[kt2][j];
   kvn1 = ltri[kt1][(i+1)%3];
   kvn2 = ltri[kt2][(j+1)%3];
   ktn1 = ltri[kt2][(j+1)%3+3];
   ktn2 = ltri[kt1][(i+1)%3+3];
   ken1 = ltri[kt2][(j+1)%3+6];
   ken2 = ltri[kt1][(i+1)%3+6];
/*
 *  If i = j = 0, the initial ltri rows are
 *
 *    kt1 = (kv1,kvn1,kvn2,kt2,ktn2,x,ke,ken2,x)
 *    kt2 = (kv2,kvn2,kvn1,kt1,ktn1,x,ke,ken1,x)
 *
 *  and the new rows, following the swap are
 *
 *    kt1 = (kv1,kvn1,kv2,ktn1,kt2,x,ken1,ke,x)
 *    kt2 = (kv2,kvn2,kv1,ktn2,kt1,x,ken2,ke,x)
 * 
 *  where the entries labeled x retain their values and need not be 
 *  referenced.
 */
   ltri[kt1][(i+2)%3] = kv2;
   ltri[kt2][(j+2)%3] = kv1;
   ltri[kt1][i+3] = ktn1;
   ltri[kt2][j+3] = ktn2;
   ltri[kt1][(i+1)%3+3] = kt2;
   ltri[kt2][(j+1)%3+3] = kt1;
   ltri[kt1][i+6] = ken1;
   ltri[kt2][j+6] = ken2;
   ltri[kt1][(i+1)%3+6] = ke;
   ltri[kt2][(j+1)%3+6] = ke;
/*
 *  Correct neighboring triangle ltri entries if necessary.
 */
   if (ktn1 >= 0) {
      if (ltri[ktn1][3] == kt2) 
         ltri[ktn1][3] = kt1;
      else if (ltri[ktn1][4] == kt2)
         ltri[ktn1][4] = kt1;
      else
         ltri[ktn1][5] = kt1;
   }
   if (ktn2 >= 0) {
      if (ltri[ktn2][3] == kt1)
         ltri[ktn2][3] = kt2;
      else if (ltri[ktn2][4] == kt1)
         ltri[ktn2][4] = kt2;
      else
         ltri[ktn2][5] = kt2;
   }
/*
 *  Adjust lvtx entries.
 */
   if (lvtx[kvn1] == kt2) lvtx[kvn1] = kt1;
   if (lvtx[kvn2] == kt1) lvtx[kvn2] = kt2;
/*
 *  Adjust ledg entries, preserving the property that the larger of 
 *  the two triangle indices is used.
 */
   ledg[ken1] = kt1;
   if (ktn1 > kt1) ledg[ken1] = ktn1;
   ledg[ken2] = kt2;
   if (ktn2 > kt2) ledg[ken2] = ktn2;
   if (!updateR) return;
/*
 *  Update lister and listtr.
 */ 
   for (k = 0; k < ner; k++) {
      if (lister[k] == ke) {
         if (!listf[kv1]  ||  !listf[kv2]) {
            ner--;
            lister[k] = lister[ner];
         }
         break;
      }
   }
   if (k == ner  &&  listf[kv1]  &&  listf[kv2]) {
      lister[ner] = ke;
      ner++;
   }
   for (k = 0; k < ntr; k++) {
      if (listtr[k] == kt1) {
         if (!listf[kv2]) {
            ntr--;
            listtr[k] = listtr[ntr];
         } 
         break;
      }
   }
   if (k == ntr  &&  listf[kv1]  &&  listf[kv2]  &&  listf[kvn1]) {
      listtr[ntr] = kt1;
      ntr++;
   }
   for (k = 0; k < ntr; k++) {
      if (listtr[k] == kt2) {
         if (!listf[kv1]) {
            ntr--;
            listtr[k] = listtr[ntr];
         }
         break;
      }
   }
   if (k == ntr  &&  listf[kv1]  &&  listf[kv2]  &&  listf[kvn2]) {
      listtr[ntr] = kt2;
      ntr++;
   }
   return;
}
/* End of swap */


/*
 *--------------------------------------------------------------------
 *
 *  swapEdges:  Mouse callback (callback triggered by a mouse button
 *              press or release with the mouse position in the 
 *              primary window) used to select an edge to be swapped.
 *
 *  This function swaps the edge selected by a left button press if
 *  the mouse position is close enough to an edge, and the edge can
 *  be swapped.
 *
 *  glmesh2 functions required:  plDistance, swap, xlate
 *
 *--------------------------------------------------------------------
 */
static void swapEdges(int button, int state, int x, int y)
{
   GLdouble wx, wy, wz;     /* Computed object coordinates */
   GLboolean intr;
   GLdouble tols;
   GLint i, ke, kt, ktn, kv0, kv1 = 0, kv2 = 0, kv3;

   if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
      xlate(x,y,0.0, &wx,&wy,&wz);
/*
 *  Compute a tolerance tols on the squared distance from the
 *  mouse position (wx,wy) to the edge.
 */
      tols = (sx <= sy) ? sx : sy;   /* tols = min(sx, sy)  */
      tols = 16.0*tols*tols;         /* 4-pixel tolerance   */
/*
 *  Loop on edges kv1-kv2.
 */
      intr = 0;
      for (kt = 0; kt < nt; kt++) {
	 for (i = 0; i < 3; i++) {
            ktn = ltri[kt][i+3];
            if (ktn > kt) continue;
	    kv1 = ltri[kt][(i+1)%3];
	    kv2 = ltri[kt][(i+2)%3];
            if ((intr = plDistance(wx, wy, vtxy[kv1][0], vtxy[kv1][1], 
                            vtxy[kv2][0], vtxy[kv2][1], tols))) break;
         }
	 if (intr) break;
      }
      if (!intr) return;
      if (ktn < 0) {
         printf("\n  A boundary edge cannot be swapped.\n");
         return;
      }
      kv0 = ltri[kt][i];
      ke = ltri[kt][i+6];
      if (ltri[ktn][6] == ke)
         kv3 = ltri[ktn][0];
      else if (ltri[ktn][7] == ke)
         kv3 = ltri[ktn][1];
      else
         kv3 = ltri[ktn][2];
/*
 *  Edge ke can be swapped iff kv1 strictly Left kv3 -> kv0 and 
 *                             kv2 strictly Left kv0 -> kv3.
 */
      if (((vtxy[kv0][0]-vtxy[kv3][0])*(vtxy[kv1][1]-vtxy[kv3][1]) >
           (vtxy[kv0][1]-vtxy[kv3][1])*(vtxy[kv1][0]-vtxy[kv3][0])       
        && (vtxy[kv3][0]-vtxy[kv0][0])*(vtxy[kv2][1]-vtxy[kv0][1]) >
           (vtxy[kv3][1]-vtxy[kv0][1])*(vtxy[kv2][0]-vtxy[kv0][0]))) {
         swap(ke,1);
         glutPostRedisplay();
      } else {
         printf("\n That edge cannot be swapped.\n");
      }
   }
   return;
}
/* End of swapEdges */


/*
 *--------------------------------------------------------------------
 *
 *  swptst:  Delaunay circumcircle test.
 *
 *  This function applies the circumcircle test to a quadrilateral
 *  defined by a pair of adjacent triangles.  The diagonal edge 
 *  (shared triangle side) should be swapped for the other diagonl if 
 *  and only if the fourth vertex is strictly interior to the circum-
 *  circle of one of the triangles (the decision is independent of the
 *  choice of triangle).  Equivalently, the diagonal is chosen to
 *  maximize the smallest of the six interior angles over the two 
 *  pairs of possible triangles (and the decision is for no swap if 
 *  the quadrilateral is not strictly convex).
 *
 *  When the four vertices are nearly cocircular (the neutral case),
 *  the preferred decision is no swap -- in order to avoid unneces-
 *  sary swaps and, more important, to avoid cycling in a sequence of
 *  sweeps (through a set of triangles) terminated by the global
 *  empty circumcircle property (no circumcircle of a triangle in
 *  the set contains a vertex in its interior).  Thus, a tolerance 
 *  swtol is used to define 'nearness' to the neutral case.
 *
 *  On input:  in1, in2, io1, and io2 are the vertex indices (indices
 *             for vtxy) of the quadrilateral with triangles defined
 *             by CCW-ordered triples (io1,io2,in1) and (io2,io1,in2)
 *             so that io1-io2 is the shared side to be replaced by
 *             edge in1-in2.
 *
 *  On output:  the return value is 1 iff the shared edge io1-io2 
 *              should be swapped for the other diagonal in1-in2.
 *              
 *  Global variable required:  vtxy
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static unsigned char swptst(int in1, int in2, int io1, int io2)
{
   double cos1, cos2, dx11, dx12, dx21, dx22, dy11, dy12, dy21, dy22,
          sin1, sin2, sin12;
   double swtol = 2.2e-15;  /* 10*(Machine precision) */
/*
 *  Compute components of the directed edges anchored at in1 and in2.
 */
   dx11 = vtxy[io1][0] - vtxy[in1][0];
   dx12 = vtxy[io2][0] - vtxy[in1][0];
   dx21 = vtxy[io1][0] - vtxy[in2][0];
   dx22 = vtxy[io2][0] - vtxy[in2][0];

   dy11 = vtxy[io1][1] - vtxy[in1][1];
   dy12 = vtxy[io2][1] - vtxy[in1][1];
   dy21 = vtxy[io1][1] - vtxy[in2][1];
   dy22 = vtxy[io2][1] - vtxy[in2][1];
/*
 *  Compute inner products, proportional cos(A1) and cos(A2), where
 *  A1 and A2 are the angles at in1 and in2, respectively.
 */
   cos1 = dx11*dx12 + dy11*dy12;
   cos2 = dx21*dx22 + dy21*dy22;
/*
 *  The diagonals should be swapped iff (A1+A2) > 180 degrees.  The 
 *  following two tests ensure numerical stability:  the decision must
 *  be FALSE when both angles are close to 0, and TRUE when both 
 *  angles are close to 180 degrees.
 */
   if (cos1 >= 0.0  &&  cos2 >= 0.0) return 0;
   if (cos1 < 0.0  &&  cos2 < 0.0) return 1;
/*
 *  Compute vector cross products (Z-components) proportional to
 *  sin(A1) and sin(A2), and sin12 proportional to sin(A1+A2).
 */
   sin1 = dx11*dy12 - dx12*dy11;
   sin2 = dx22*dy21 - dx21*dy22;
   sin12 = sin1*cos2 + cos1*sin2;
   return (sin12 < -swtol);
}
/* End of swptst */


/*
 *--------------------------------------------------------------------
 *
 *  trmtst:  Triangulation data structure integrity test.
 *
 *  This function tests for the following properties:
 *
 *    1)  Counts nc, nb, nv, nn, nt, and ne are consistent.
 *    2)  Vertex indices in ltri, lcrv, and lnrv are in the range 0 
 *        to nv-1.
 *    3)  Edge indices in ltri are in the range 0 to ne-1.
 *    4)  Triangle indices in ltri, lvtx, and ledg are in the range
 *        0 to nt-1.
 *    5)  Triangle ledg[ke] contains edge ke, and ledg[ke] is the 
 *        larger of the two relevant triangle indices for all ke.
 *    6)  Triangle lvtx[kv] contains vertex kv for all kv.
 *    7)  Adjacency information is consistent.
 *
 *  On input:
 *
 *      nc,nb,nv,nn,nt,ne = Numbers of boundary curves, boundary 
 *                          edges, vertices, non-removable vertices,
 *                          triangles, and edges, respectively.
 *
 *      ltri,lvtx,ledg,lnrv,lcrv = Triangle list, vertex list, edge
 *                                 list, non-removable vertex list,
 *                                 and boundary curve list.
 *
 *  If an integrity check is not satisfied, an error message is
 *  printed, and the program is terminated.  Otherwise a message
 *  indicating success is printed.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void trmtst(int nc, int nb, int nv, int nn, int nt, int ne, 
                   int ltri[][9], int lvtx[], int ledg[], int lnrv[],
                   int lcrv[])
{
   int i, j, ke, kt, kt0, kv;
/*
 *  Test for consistency of counts.
 */
   if (nb < 3*nc  ||  nn < 0  ||  nn > nv  ||
       ne != 3*nv-nb+3*nc-6  ||
       nt != 2*nv-nb+2*nc-4) {
      printf("glmesh2 (trmtst):  Inconsistent counts:  nv = %d, "
             "nb = %d, nc = %d, nt = %d, ne = %d\n", nv, nb, nc, nt, 
             ne);
      exit(1);
   }
/*
 *  Test ltri values.
 */
   for (kt = 0; kt < nt; kt++) {
      if (ltri[kt][0] < 0  ||  ltri[kt][0] >= nv  ||
          ltri[kt][1] < 0  ||  ltri[kt][1] >= nv  ||
          ltri[kt][2] < 0  ||  ltri[kt][2] >= nv  ||
          ltri[kt][3] < -1  ||  ltri[kt][3] >= nt  ||
          ltri[kt][4] < -1  ||  ltri[kt][4] >= nt  ||
          ltri[kt][5] < -1  ||  ltri[kt][5] >= nt  ||
          ltri[kt][6] < 0  ||  ltri[kt][6] >= ne  ||
          ltri[kt][7] < 0  ||  ltri[kt][7] >= ne  ||
          ltri[kt][8] < 0  ||  ltri[kt][8] >= ne) {
         printf("glmesh2 (trmtst):  Invalid ltri[kt] entry:  kt = "
                "%d\n", kt);
         exit(1);
      }
   }
/*
 *  Test ledg values.
 */
   for (ke = 0; ke < ne; ke++) {
      kt = ledg[ke];
      if (kt < 0  ||  kt >= nt) {
         printf("glmesh2 (trmtst):  Invalid kt = ledg[ke] entry:  "
                "ke = %d, kt = %d\n", ke, kt);
         exit(1);
      }
      if (ltri[kt][6] == ke)
         i = 3;
      else if (ltri[kt][7] == ke)
         i = 4;
      else if (ltri[kt][8] == ke)
         i = 5;
      else {
         printf("glmesh2 (trmtst):  Invalid kt = ledg[ke] or "
                "ltri[kt] entry:  ke = %d, kt = %d\n", ke, kt);
         exit(1);
      }
      kt0 = kt;
      kt = ltri[kt][i];
      if (kt >= kt0) {
         printf("glmesh2 (trmtst):  Invalid kt = ledg[ke] entry (not"
                " the smaller triangle index:  ke = %d, kt = %d\n", 
                ke, kt0);
         exit(1);
      }
      if (kt < 0) continue;
      if (ltri[kt][6] == ke)
         i = 3;
      else if (ltri[kt][7] == ke)
         i = 4;
      else if (ltri[kt][8] == ke)
         i = 5;
      else {
         printf("glmesh2 (trmtst):  Invalid ledg[ke] or ltri[kt] "
                "entry:  ke = %d, kt = %d\n", ke, kt);
         exit(1);
      }
      if (ltri[kt][i] != kt0) {
         printf("glmesh2 (trmtst):  Inconsistency between ltri[kt1] "
                "and ltri[kt2], kt1 = %d, kt2 = %d\n", kt0, kt);
         exit(1);
      }
   }
/*
 *  Test lvtx values.
 */
   for (kv = 0; kv < nv; kv++) {
      kt = lvtx[kv];
      if (kt < 0  ||  kt >= nt) {
         printf("glmesh2 (trmtst):  Invalid kt = lvtx[kv] entry:  "
                "kv = %d, kt = %d\n", kv, kt);
         exit(1);
      }
      if (ltri[kt][0] == kv)
         i = 0;
      else if (ltri[kt][1] == kv) 
         i = 1;
      else if (ltri[kt][2] == kv)
         i = 2;
      else {
         printf("glmesh2 (trmtst):  Invalid kt = lvtx[kv] or "
                "ltri[kt] entry:  kv = %d, kt = %d\n", kv, kt);
         exit(1);
      }
   }
/*
 *  Test lnrv values.
 */
   for (i = 0; i < nn; i++) {
      kv = lnrv[i];
      if (kv < 0  ||  kv >= nv) {
         printf("glmesh2 (trmtst):  Invalid kv = lnrv[k] entry:  k ="
                " %d, kv = %d\n", i, kv);
         exit(1);
      }
      kt = lvtx[kv];
      if (ltri[kt][0] == kv)
         j = 2;
      else if (ltri[kt][1] == kv)
         j = 0;
      else
         j = 1;
      if (ltri[kt][j+3] >= 0) {
         printf("glmesh2 (trmtst):  Vertex kv = lnrv[k] is not a "
                "boundary vertex:  k = %d, kv = %d\n", i, kv);
         exit(1);
      }
   }
/*
 *  Test lcrv values.
 */
   for (i = 0; i < nc; i++) {
      kv = lcrv[i];
      if (kv < 0  ||  kv >= nv) {
         printf("glmesh2 (trmtst):  Invalid kv = lcrv[k] entry:  k ="
                " %d, kv = %d\n", i, kv);
         exit(1);
      }
      kt = lvtx[kv];
      if (ltri[kt][0] == kv)
         j = 2;
      else if (ltri[kt][1] == kv)
         j = 0;
      else
         j = 1;
      if (ltri[kt][j+3] >= 0) {
         printf("glmesh2 (trmtst):  Vertex kv = lcrv[k] is not a "
                "boundary vertex:  k = %d, kv = %d\n", i, kv);
         exit(1);
      }
   }
/*
 *  No error encountered.
 */
   printf("\n  Triangulation integrity tests passed.\n");
   return;
}
/* End of trmtst */


/*
 *--------------------------------------------------------------------
 *
 *  trprint:  Print triangle mesh data structure on standard output
 *            unit.  Triangle shape quality statistics, and the lists
 *            associated with R are also printed.
 *
 *  On input:
 *
 *      nc,nb,nv,nn,nt,ne = Numbers of boundary curves, boundary 
 *                          edges, vertices, non-removable vertices,
 *                          triangles, and edges, respectively.
 *
 *      vtxy = Type double array dimensioned nv by 2 containing the
 *             Cartesian coordinates of the vertices.
 *
 *      uv = Type double array dimensioned nv by 2 containing pairs
 *           of function values (u,v) at the vertices.
 *
 *      ltri,lvtx,ledg,lnrv,lcrv = Triangle list, vertex list, edge 
 *                                 list, non-removable vertex list,
 *                                 and boundary curve list.
 *
 *  glmesh2 functions called:  initR, trqual1, trqual2
 *
 *--------------------------------------------------------------------
 */
static void trprint(int nc, int nb, int nv, int nn, int nt, int ne, 
                    double vtxy[][2], double uv[][2], int ltri[][9], 
                    int lvtx[], int ledg[], int lnrv[], int lcrv[])
{
   char *s;
   char *separator = "\n________________________________________"
                     "________________________________________\n";
   double pfair, pgood, q, qmax, qmean, qmin;
   int i, kt, kv, kv1, kv2, kv3;
/*
 *  Store lists associated with R if necessary.
 */
   if (newR) initR();
/*
 *  Print title and counts.
 */
   printf("%s", separator);
   printf("         Triangle Mesh Data Structure\n\n"
          "  Nc = %d boundary curves\n"
          "  Nb = %d boundary edges (boundary vertices)\n"
          "  Nv = %d vertices\n"
          "  Nn = %d non-removable boundary vertices\n"
          "  Nt = %d triangles\n"
          "  Ne = %d edges\n\n", nc, nb, nv, nn, nt, ne);
/*
 *  Print the vertex coordinates, (u,v) pairs, vertex list lvtx,
 *  and *'s after the non-removable vertices.
 */
   printf("                                  Vertex List\n\n"
          "  Vertex           x           y           u           v"
          "    lvtx  Non-removable\n\n");
   for (kv = 0; kv < nv; kv++) {
      s = "";
      for (i = 0; i < nn; i++) {
         if (lnrv[i] == kv) {
            s = "*";
            break;
         }
      }
      printf("  %6.1d  %10.3e  %10.3e  %10.3e  %10.3e  %6.1d"
             "          %s\n", kv, vtxy[kv][0], vtxy[kv][1], 
             uv[kv][0], uv[kv][1], lvtx[kv], s); 
   }
/*
 *  Print the triangle list.
 */
   printf("\n\n                              Triangle List\n\n"
          "      Vertex Indices       Neighboring Triangles      "
          "Edge Indices\n\n");
   for (i = 0; i < nt; i++)
      printf("  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d"
             "  %6.1d  %6.1d\n", ltri[i][0], ltri[i][1], ltri[i][2],
             ltri[i][3], ltri[i][4], ltri[i][5], ltri[i][6], 
             ltri[i][7], ltri[i][8]);
/*
 *  Print the edge list.
 */
   printf("\n\n  Edge list (triangle index for each edge):\n\n");
   for (i = 0; i < ne-8; i += 9) 
      printf("  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d"
             "  %6.1d  %6.1d\n", ledg[i], ledg[i+1], ledg[i+2],
             ledg[i+3], ledg[i+4], ledg[i+5], ledg[i+6], ledg[i+7],
             ledg[i+8]);
   for (i = ne-(ne%9); i < ne; i++) 
      printf("  %6.1d", ledg[i]);
   printf("\n");
/*
 *  Print the boundary curve list.
 */
   printf("\n\n  Boundary curve list (boundary vertex index for each"
          " curve):\n\n");
   for (i = 0; i < nc-8; i += 9) 
      printf("  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d  %6.1d"
             "  %6.1d  %6.1d\n", lcrv[i], lcrv[i+1], lcrv[i+2],
             lcrv[i+3], lcrv[i+4], lcrv[i+5], lcrv[i+6], lcrv[i+7],
             lcrv[i+8]);
   for (i = nc-(nc%9); i < nc; i++) 
      printf("  %6.1d", lcrv[i]);
   printf("\n");
/*
 *  Compute and print shape-quality statistics.
 */
   qmin = 2.0;
   qmean = 0.0;
   qmax = 0.0;
   pfair = 0.0;
   pgood = 0.0;
   for (kt = 0; kt < nt; kt++) {
      kv1 = ltri[kt][0];
      kv2 = ltri[kt][1];
      kv3 = ltri[kt][2];
      q = trqual1(vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0],
                  vtxy[kv2][1], vtxy[kv3][0], vtxy[kv3][1]);
      if (q < qmin) qmin = q;
      qmean += q;
      if (q > qmax) qmax = q;
      if (q > 0.6) pfair += 1.0;
      if (q > 0.866) pgood += 1.0;
   }
   qmean /= (GLdouble) nt;
   pfair = 100.0*pfair/(GLdouble) nt;
   pgood = 100.0*pgood/(GLdouble) nt;
   printf("\n  Triangle shape quality q1 = "
          "Area/Mean-squared-side-length (normalized):\n");
   printf("     q1_min = %5.3f      Percent poor (q1 < 0.6) = %6.2f\n",
          qmin, 100.0-pfair);
   printf("     q1_mean = %5.3f     Percent fair (q1 >= 0.6) = %6.2f\n",
          qmean, pfair);
   printf("     q1_max = %5.3f      Percent good (q1 >= 0.866) = %6.2f\n",
          qmax, pgood);

   qmin = 2.0;
   qmean = 0.0;
   qmax = 0.0;
   pfair = 0.0;
   pgood = 0.0;
   for (kt = 0; kt < nt; kt++) {
      kv1 = ltri[kt][0];
      kv2 = ltri[kt][1];
      kv3 = ltri[kt][2];
      q = trqual2(vtxy[kv1][0], vtxy[kv1][1], vtxy[kv2][0],
                  vtxy[kv2][1], vtxy[kv3][0], vtxy[kv3][1]);
      if (q < qmin) qmin = q;
      qmean += q;
      if (q > qmax) qmax = q;
      if (q > 0.5) pfair += 1.0;
      if (q > 0.7) pgood += 1.0;
   }
   qmean /= (GLdouble) nt;
   pfair = 100.0*pfair/(GLdouble) nt;
   pgood = 100.0*pgood/(GLdouble) nt;
   printf("\n  Triangle shape quality q2 = "
          "Inradius/Circumradius (normalized):\n");
   printf("     q2_min = %5.3f      Percent poor (q2 < 0.5) = %6.2f\n",
          qmin, 100.0-pfair);
   printf("     q2_mean = %5.3f     Percent fair (q2 >= 0.5) = %6.2f\n",
          qmean, pfair);
   printf("     q2_max = %5.3f      Percent good (q2 >= 0.7) = %6.2f\n",
          qmax, pgood);
/*
 *  Print the lists associated with R.
 */
   printf("\n\n  Number of triangles in R, ntr = %d, listtr =\n", ntr);
   for (i = 0; i < ntr; i++) printf("%6.1d  ", listtr[i]);
   printf("\n\n  Number of edges in R, ner = %d, lister =\n", ner);
   for (i = 0; i < ner; i++) printf("%6.1d  ", lister[i]);
   printf("\n\n  Number of vertices in R, nvr = %d, listvr =\n", nvr);
   for (i = 0; i < nvr; i++) printf("%6.1d  ", listvr[i]);
   printf("\n\n  Vertex flags, listf =\n  ");
   for (i = 0; i < nv; i++) printf("%u  ", listf[i]);
   printf("%s", separator);
   fflush(stdout);
   return;
}
/* End of trprint */


/*
 *--------------------------------------------------------------------
 *
 *  trqual1:  Compute a measure of triangle quality,
 *
 *                 q = 4*sqrt(3)*A/(s1*s1 + s2*s2 + s3*s3),
 *
 *            where A is the triangle area, and s1, s2, and s3 are
 *            the side lengths.  The scaling is chosen so that an
 *            equilateral triangle has q = 1, which is optimal.  
 *
 *  We define three categories:
 *
 *  1)  Good:  q >= sqrt(3)/2.  The largest angle is at most pi/2, and
 *             the smallest angle is at least arccos(.8) in this case.
 *
 *  2)  Fair:  0.6 <= q < sqrt(3)/2.  The largest angle is at most
 *             2*pi/3, and the smallest angle is at least 
 *             arccos(13/14) in this case.
 *
 *  3)  Poor:  q < 0.6.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static double trqual1(double x1, double y1, double x2, double y2,
		      double x3, double y3)
{
   double a, dx1, dx2, dx3, dy1, dy2, dy3, s1, s2, s3;
/*
 *  Compute squared side lengths s1, s2, and s3.
 */   
   dx1 = x3-x2;
   dy1 = y3-y2;
   dx2 = x1-x3;
   dy2 = y1-y3;
   dx3 = x2-x1;
   dy3 = y2-y1;
   s1 = dx1*dx1 + dy1*dy1;
   s2 = dx2*dx2 + dy2*dy2;
   s3 = dx3*dx3 + dy3*dy3;
/*
 *  Compute a = 2*area.
 */
   a = fabs( dx3*dy2 - dx2*dy3 );
   return (2.0*sqrt(3.0)*a/(s1 + s2 + s3));
}	
/* End of trqual1 */


/*
 *--------------------------------------------------------------------
 *
 *  trqual2:  Compute a measure of triangle quality,
 *
 *              q = 2*r/R = (s2+s3-s1)*(s3+s1-s2)*(s1+s2-s3)/s1*s2*s3,
 *
 *            where r and R denote the inradius and circumradius of 
 *            the triangle, and s1, s2, and s3 are the side lengths.  
 *            The value ranges from 0 for a degenerate triangle to 1
 *            for an equilateral triangle.
 *
 *  We define three categories:
 *
 *  1)  Good:  q >= 0.7
 *
 *  2)  Fair:  0.5 <= q < 0.7.
 *
 *  3)  Poor:  q < 0.5.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static double trqual2(double x1, double y1, double x2, double y2,
		      double x3, double y3)
{
   double dx1, dx2, dx3, dy1, dy2, dy3, s1, s2, s3;
/*
 *  Compute side lengths s1, s2, and s3.
 */   
   dx1 = x3-x2;
   dy1 = y3-y2;
   dx2 = x1-x3;
   dy2 = y1-y3;
   dx3 = x2-x1;
   dy3 = y2-y1;
   s1 = sqrt(dx1*dx1 + dy1*dy1);
   s2 = sqrt(dx2*dx2 + dy2*dy2);
   s3 = sqrt(dx3*dx3 + dy3*dy3);
 
   return ((s2+s3-s1)*(s3+s1-s2)*(s1+s2-s3)/(s1*s2*s3));
}	
/* End of trqual2 */


/*
 *--------------------------------------------------------------------
 *
 *  trwrite:  Output the triangle mesh to a file.
 *
 *  On input:
 *
 *      fname = File name for output.
 *
 *      nc,nb,nv,nt,ne = Numbers of boundary curves, boundary edges,
 *                       vertices, triangles, and edges, respectively.
 *
 *      vtxy = Type double array dimensioned nv by 2 containing the
 *             Cartesian coordinates of the vertices.
 *
 *      uv = Type double array dimensioned nv by 2 containing pairs
 *           of function values (u,v) at the vertices.
 *
 *      ltri = Triangle list.
 *
 *  glmesh2 functions called:  None
 *
 *--------------------------------------------------------------------
 */
static void trwrite(char *fname, int nc, int nb, int nv, int nt, 
                    int ne, double vtxy[][2], double uv[][2],
                    int ltri[][9])
{
   FILE *fpw;             /* File pointer */
   time_t curtime;        /* Current calendar time */
   struct tm *loctime;    /* Local time representation */
   int  linsiz = 160;     /* Buffer size */
   char line[160];        /* Buffer for output line */
   int kt, kv, nbc = 0;
/*
 *  Open file fname with write access.
 */
   fpw = fopen(fname, "w");
   if (fpw == NULL) {
      printf("glmesh2:  Cannot open file %s.\n", fname);
      exit(1);
   }
/*
 *  Output header comments.
 */
   curtime = time (NULL);
   loctime = localtime (&curtime);
   strftime(line,linsiz,"#\n# Data file created by glmesh2 on "
            "%b %d, %Y\n", loctime);
   fputs(line,fpw);
   fprintf(fpw,"# Nc = %d, Nb = %d, Nv = %d, Nt = %d, Ne = %d\n",
            nc, nb, nv, nt, ne);
   fprintf(fpw,"#\n# Vertex list:  Nv followed by Nv (x,y,u,v) "
               "values\n");
/*
 *  Output vertex list.
 */
   fprintf(fpw," %d\n", nv);
   for (kv = 0; kv < nv; kv++) {
      fprintf(fpw," %19.12e %19.12e %19.12e %19.12e\n", vtxy[kv][0], 
              vtxy[kv][1], uv[kv][0], uv[kv][1]);
   }
/*
 *  Output the triangle list.
 */
   fprintf(fpw,"#\n# Triangle list:  Nt followed by Nt CCW-ordered"
               " triples of vertex indices\n");
   fprintf(fpw," %d\n", nt);
   for (kt = 0; kt < nt; kt++) {
      fprintf(fpw," %9d %9d %9d\n", ltri[kt][0], ltri[kt][1],
              ltri[kt][2]);
   }
/*
 *  Output a boundary segment list.
 */
   fprintf(fpw,"#\n# Boundary segment list:  Nbc followed by Nbc"
               " triples (type,ke1,ke2)\n");
   fprintf(fpw," %d\n", nbc);
   return;
}
/* End of trwrite */


/*
 *--------------------------------------------------------------------
 *
 *  unrefine:  Unrefine the mesh by reversing the previous refinement
 *             by refineRe or refineRt unless there was no previous
 *             refinement or the previous refinement was rendered non-
 *             reversible by deletion of a vertex (in Function 
 *             coarsenRv or removeRt) or by a change in lnrv.
 *
 *  glmesh2 function called:  delvtx
 *
 *--------------------------------------------------------------------
 */
static void unrefine(void)
{
   int kv, nfail, nv0;

   if (nref <= 0) {
      printf("\n  Previous refinement cannot be reversed.\n");
      return;
   }
/*
 *  Pop rstack into nv0.
 */
   nref--;
   nv0 = rstack[nref];
/*
 *  Delete vertex sequence nv-1, nv-2, ..., nv0.
 */
   nfail = 0;
   for (kv = nv-1; kv >= nv0; kv--) { 
      if (!delvtx(kv)) nfail++;
   }
   return;
}
/* End of unrefine */


/*
 *--------------------------------------------------------------------
 *
 *  xlate:  Translate glut window coordinates (wx,wy), with depth wz, 
 *          to object coordinates (x,y,z) by inverting the mappings 
 *          associated with the modelview, projection, and viewport 
 *          matrices.
 *
 *  On input:  wx and wy are the integer window-relative coordinates
 *             in the range 0 to w-1 and 0 to h-1, respectively, for
 *             window width w and height h, and wz is the depth in
 *             the range 0 to 1.0.  The origin is at the upper left
 *             corner of the window.
 *
 *  On output:  x, y, and z are the objects coordinates.
 *
 *  Global variables required:  None
 *
 *  glmesh2 functions required:  None
 *
 *--------------------------------------------------------------------
 */
static void xlate(GLint wx, GLint wy, GLdouble wz, 
                  GLdouble *x, GLdouble *y, GLdouble *z)
{
   GLdouble mvmatrix[16];    /* Modelview matrix */
   GLdouble projmatrix[16];  /* Projection matrix */
   GLint viewport[4];        /* Components of Viewport mapping */
   GLint yr;                 /* OpenGL y coordinate */

   glGetIntegerv(GL_VIEWPORT, viewport);
   glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
   glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);
/*
 *  Compute yr = h-1-wy for window height h.
 */
   yr = glutGet(GLUT_WINDOW_HEIGHT) - 1 - wy;
   gluUnProject((GLdouble) wx, (GLdouble) yr, wz, mvmatrix, 
                projmatrix, viewport, x, y, z);
   return;
}
/* End of xlate */


/*
 *--------------------------------------------------------------------
 *
 *  main:  Input data, create windows, create a menu, register call-
 *         back functions, and enter a loop.
 *
 *  glmesh2 functions called:  getCorners, init1, inputData, makeMenu,
 *                             storeV
 *
 *  glmesh2 callback functions registered:  display, dlgDisplay, 
 *                   dlgkey, dlgMouse, dlgReshape, key, motion, mouse,
 *                   passiveMotion, reshape, specialKey
 *
 *--------------------------------------------------------------------
 */
int main(int argc, char *argv[])
{
   inputData(argc, argv);        /* Input data */
   storeV();                     /* Initialize global parameters */
   getCorners();                 /* Initialize lnrv to corners */
   glutInit(&argc, argv);        /* Initialize Glut library */
/*
 *  Create the primary window (surface display):  window 1.
 */
   glutInitWindowPosition(50, 50);
   glutInitWindowSize(500, 500);    /* Size in pixels */
   glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
   window[0] = glutCreateWindow(argv[1]);      /* Create window 1 */

   makeMenu();                   /* Create a menu attached
                                    to the right button */

   glutKeyboardFunc(key);        /* Register callbacks */
   glutSpecialFunc(specialKey);
   glutMouseFunc(mouse);
   glutMotionFunc(motion);
   glutPassiveMotionFunc(passiveMotion);
   glutReshapeFunc(reshape); 
   glutDisplayFunc(display);

   init1();                      /* One-time initialization */
/*
 *  Create the secondary window (dialog):  window 2.
 */
   glutInitWindowPosition(500, 250);
   glutInitWindowSize(400, 200);         /* Size in pixels */
   glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE);
   window[1] = glutCreateWindow("Dialog");  /* Create window 2 */
   glClearColor (0.5, 0.5, 0.5, 1.0);       /* Grey background */

   glutKeyboardFunc(dlgKey);       /* Register callbacks */
   glutMouseFunc(dlgMouse);
   glutReshapeFunc(dlgReshape);
   glutDisplayFunc(dlgDisplay);
   glutHideWindow();             /* Initially hidden */

   glutSetWindow(window[0]);     /* Make window 1 current */
   glutMainLoop();
   return 0;
}
/* End of main */
